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ABSTRACT 

The  accuracy  of  linear  multistep  formulas  suitable  for  stiff 
differential   systems   is  limited.     Greater  accuracy  can  be  attained  by 
including  higher  derivatives   in  the  formula,  but  this  is  not  practical 
for  all   problems.     However  it  is  possible  to  duplicate  the  absolute 
stability  region  for  any  given  m-derivative  multistep  formula  by  taking 
a  linear  combination  of  m  multistep  formulas.     These  blended  formulas 
are  similar  to  Lambert  and  Sigurdsson's  linear  multistep  formulas  with 
variable  matrix  coefficients,  but  the  approach  is  different.     Imple- 
mentation details  and  numerical    results  are  presented  for  a  variable- 
order,   variable-step  blend  of  the  Adams-Moulton  and  the  backward 
differentiation  formulas. 


910. 

1.      INTRODUCTION 

Currently  the  most  popular  codes  for  solving  systems  of  stiff 
ordinary  differential   equations   are  based  on  the  backward  differentiation 
formulas,  for  example  Gear  [9],  Bray  ton  et  al    [  1  ],  Hindmarsh  [12],  and 
Byrne  and  Hindmarsh  [3].     The  biggest  drawback  [7,  p.   33]  of  these  formulas 
is   the  poor  stability  properties  of  the  higher  order  formulas  when  the 
Jacobian  matrix  for  the  problem  has  eigenvalues  close  to  the  imaginary 
axis.     This  difficulty  can  be  overcome  by  restricting  the  order  of  the 
formulas;  however  these  lower  order  formulas  are  often  fairly  inefficient 
because  of  their  limited  accuracy. 

There  are  conceivably  three  situations  in  which  more  accurate 
stiffly  stable  formulas  would  be  helpful: 

(i)   for  nonstiff  problems.     The  casual    user  of  a  program  from  a 

subroutine  library  should  not  be  reauired  to  know  whether  or  not 
his   problem  is  stiff.     And  there  is  presently  no  auick  and 
reliable  way  for  automatically  determining  whether  or  not  a 
problem  is  stiff.     Therefore  it  may  be  reasonable  to  provide  a 
stiff  formula  as  the  default.      (If  efficiency  were  very  imnortant, 
the  user  should  not  be  considered  as   "casual".) 
(ii)   for  stiff  problems  with  moderate  to  high  accuracy  requirements, 
(iii)   for  stiff  problems  whose  Jacobian  matrix  has  large  eigenvalues 
near  the  imaginary  axis.     These  are  the  problems   for  which   the 
backward  differentiation  formulas   are  often  inadeouate. 

One  can  get  more  accurate  stiffly  stable  formulas  for  a  system 

&-f(t.y> 


of  s  differential   equations  by  including  derivatives  of  f  in  the  formula: 

k  m       .     k  f.n 

z     o.  yM   .  =     z     hJ     z     b  ••   f     .    ' 
1=0     n     n_1       J-l         1-0     J1     n"n 

where  f^       =  f  ^n-1 '  yn-i  ^     Here  f        denotes   the  total    derivative 

of  f  defined  recursively  by 

f<°>  =  f     , 

fU+l)    =    9_  f(£)    +    f   3_  fU) 

at  3y 

Formulas  of  this  type  have  been  proposed  by  Enright  [5]  and  Rrown   [2]. 
The  presence  of  higher  derivatives  makes  it  necessary  to  reauire   the  user 
to  provide   routines   for  their  evaluation  or  better  still   to  have  these 
routines  generated  symbolically  from  f(t,  y).     Unfortunately  there  are 
differential   equations   for  which  symbolic  differentiation  is  not  possible, 
and  in  these  cases   the  user  may  be  unwilling  (see   [18])   or  unable  to 
provide  analytic  Jacobians.     Thus  the  biggest  drawback  of  these 
multi derivative  formulas  is   their  limited  applicability. 

Stiff  problems  require  implicit  formulas,  and  hence  systems  of 
nonlinear  equations  must  be  solved.     This  is  usually  accomplished  by  some 
kind  of  modified  Newton  iteration,  which  means  that  an  approximation  to 
the  Jacobian  is  available  during  the  integration  of  a  stiff  system.     It  may 
be  advantageous   to  make  further  use  of  the  Jacobian  in  our  method.     We 
cannot  use  it  to  improve  the  order  of  accuracy  of  the  method  because  it  may 
not  be  exact,  but  we  can  use  it  to  enhance  the  stability  properties.     This 
idea  has  been  pursued  by  Lambert  [15]  and  by  Lambert  and  Sigurdsson   [16], 
who   introduce  the   class  of  linear  multistep  formulas  with  variable  matrix 
coefficients.     The  formulas  presented  in   this  paper  are  of  this   type,   but 
the  approach  is   different. 


The  basic  idea  is   to  take  a  linear  mul  tiderivati  ve  multistep  formula 
and  convert  it  into  a  variable  coefficient  linear  multistep  formula  so  that 
the  region  of  absolute  stability  and  the  order  of  accuracy  is  unaffected. 
There  is  no   reason  to  suspect  that  this   transformation  would  either  improve 
or  degrade  the  performance  of  the  formula  in  the  case  of  smooth  nroblems, 
although  it  should  be  beneficial    for  problems  with  nonexistent  higher 
derivatives . 

As   an  example,   the   (k  +   1  )st  order  Enri  ght 'formula   [5]   is 
transformed  into  a  blend  of  the   (k  +  l)st  order  Adams-Moul  ton  formula  and 
the  k-th  order  backward  differentiation  formula: 

{EnrightF(k+1)}  ->  {Af1F(k+l)}  +  k  y*  h  J{BDF(k)}. 

Mere  J  is   a  Jacobian  approximation  and  y*  is   a  negative   constant  defined 

K 

in  [11,  d.    195]  and  [10,  o.   11 3].     Py  adjusting   the  coefficient  k  y*,  it 

is  possible  to  widen   the  stability  region  without  sacrificing  accuracv. 

f 
Blended     formulas  modified  in   this  way  have  been  put  into   the  CACM  alaorithm 

DIFSUB   (Gear  [9])   and  compared  with   the  BDF's    (see   §6).      The  implementation 
of  these  formulas   required  the   derivation  of  the  flordsieck  representation 
(§3)   and  the  development  of  local   error  estimators   (§5). 

There  is   a  difficulty  associated  with  thp  use  of  either 
mul  ti  deri  vati  ve  or  variable  matrix  coefficient  multistep  formulas:     the 
use  of  modified  Newton  iterations   to  solve  the  implicit  equations   requires 
"inverting"  polynomials  in  the  Jacobian  matrix.     This  would  be  undesirable 
if  sparse  matrix  techniques   are  being  used,  because  matrix  mul  tinl  i  cations 
can  considerably  reduce  sparseness.     Two  ways  of  avoiding  matrix  multipli- 
cations have  been  suggested  in  the  literature  for  the  case  m  =  ?.,     "illoughby 
t  A  term  suggested  by  C.  W.   Gear 


[6,   p.   102]  suggests   doing  the  LU  decomposition  in  conplex  arithmetic. 
Enright  [6]  suggests   choosing  the  parameters  of  the  formula  so  that  the 
quadratic  in  h  J  becomes   a  perfect  square;   the   computation   reouired  is 
one  LU  decomposition   (in  real   arithmetic)  and  Uio  backsolvos.     A  third 
alternative  is  suggested  in   §4,  which  also   requires  one  decomoosi  tion 
and  two  backsolves.     These  last  two  methods  of  avoiding  matrix  multipli- 
cations  also  avoid  the  multiplication  of  a  vector  by  a  matrix,  which  is 

(k) 
suggested  by  the  appearance  of  J{BDFVIW}  in  the  AMF-BDF  blend. 


2.      DERIVATION  OF  BLENDED  FORMULAS 

We  discuss  two  ways  of  converting  a  linear  n-deri  vative  multisteo 
formula  into  a  linear  combination  of  m  linear  multisteo  formulas: 
(i)  mixed  order  approach) 
(ii)   pure  order  approach. 
For  each  approach  an  example  is  presented  followed  by  a  brief  discussion 
of  the  general   case. 

Consider  Enright's  third  order  formula: 

2h    ,        h  f  h2  fl 

y    -y    -i=—T"f   +o-f    -i  -  r—  f    • 

7n      ■'n-l         3     n       3     n-1       6       n 


The  associated  linear  operator, 


,    2h       ',^    .    h  ..,/.       ^       h2 


Ln  y(t)  =  -y(t)  +  y(t  -  h)  +  ^y'(t)  +  ^y'(t  -  h)   -  £-y"(t)    , 

satisfies  L,    y(t)  =  0(h   ).     The  main  idea  is   to  snlit  L.    into  a  sum 
L.    +  h  -rr  Lr  where  L,    and  L,'  are  1-derivative  operators   (that  is,  the 
second  derivative  of  y(t)   is  not  present).     Me  choose  L,    so   that 
lI  y(t)  =  0(h4),  which  then  forces  L2  y(t)  =  0(h3).     Th^  conditions  on 
L,    imply  that 

lJ  y(t)  =  -y(t)  +  y(t  -  h)  +  h  30  y'(t)  +  h  b]  y'(t  -  h) 
+  ...  +  h  ek  y*(t  -  k  h). 

where  the  e's   are  chosen  so  that  L    v(t)   =  0(h   ).     There  is   a  unique 
choice  which  minimizes  the  stepn umber  k: 

Ljy(t)  =  -y(t)  +  y(t  -  h)  +  f|y'(t)  +  ^|y'(t  -  h)  -  ^-v'(t  -  2h), 

which  corresponds  to  the  third  order  Adams-foul  ton  formula.     This   choice 


of  L,    forces 
h 

Lj;  y(t)  =  -  g-{-  |y(t)  +  2y(t  -  h)  -  l-y(t  -  2h)  +  h  y'(t)}  . 

The  expression  in  braces   corresponds  to  the  second  order  backward  differen- 

tiation  formula.     Having  obtained  our  decomposition  L.    +  h  tt  L.    ,  we 

1  2 

construct  the  blended  formula  by  replacing  L,    and  L,    with  their  corresponding 

formulas  and  by  replacing  -rr  wi  th   the  Jacobian  matrix  J  =   fv(t«,  y«): 

In  an  actual    implementation  J  would  be  updated  whenever  f    was  reevaluated. 

Two  important  properties  of  this  blended  formula  should  be  noted. 
First,  the  blended  formula  is  identical  to  Enright's  formula  for  the  test 
equation  y'   =  Ay,   and  hence  the  region  of  absolute  stability 

R  =  {hx:     y     ->  0  as  n  •>  «  for  y*   =  Ay} 
is  the  same  for  both  formulas.     Second,   the  normalized  truncation  error  is 

for  some  t  and  x  between  t  -  2h   and  t,   and  so  the  formula  is  of  order  three. 
Moreover,   the  terms  in  braces   tend  to  cancel   each  other,   particularly  for 
the  equation  y1   =  Jy  for  which  the  truncation  error  would  be 

[1   -  hJ]-l{^!y(4)(t)  +  0(h5)}. 

Here  1    denotes   the  s   x  s  identity  matrix.     This   convention  of  using  scalars 
to  represent  scalar  multiples  of  the  s   *  s  identity  matrix  is   used  throughout 
the  paper. 

Having  presented  an  example,  let  us  briefly  look  at  the  general 
case  where 


k  m       .     k  f.v 

L.    y(t)  =  -     i     a.  y(t  -  ih)  +     z     hJ     i     g..  yu,{t  -  ih)   . 
n  i=0     1  j=l         i=0     J1 

Assume  the  order  p  5  m  -  2,  and  let  q  =  max{k,  p  -  1}.     We  write 
1      .  .1   .  l     d  ,Z  ,  ,   ,  m-1     d        .  m 

Lh  -  Lh +  h  dt  Lh +  •••  +  h     ^n"Lh 

where 

Ljj  y(t)   =    -a3Q  y(t)   -    ...    -  aj|  y(t  -  qh)   +  h   Bfj  y'(t)   +   ...   +  h   s;J  y'(t  -  qh) 

"I  O  ml 

This   decomposition  is   uniquely  specified  by  choosing  L,  L,  ,    ...,  L  ~ 


(in  that  order)  so  that 


and 


ii  =  0        ,        p  -  j   <   i    $  q 


il  y(t)  =  o(hp+2~j)  . 

10  m 

Having  obtained  L,  L,    ...,  L.  ,  we  define  the  mixed  order  blended  formula 
to  be 

m  .  ,     q  .  '      . 

Z     (hJ)J_l     E     (-a]  y         +  h  ^   f     •)   =  0   .  (2.1) 

j=l  i=0         n     n_1  1     n_1 

In  the  second  approach  all   the  formulas  in  the  blend  are  of  the 

same  order.     Again  we  illustrate  by  means  of  Enright's   third  order  formula. 

1  A       0  1  O 

The  operator  L.    is  split  into  a  sum  L,    +  h  -rr  L    where  L,    and  L^  are 
1 -derivative  operators.     Except  this  time  the  decomposition   is  uniquely 
specified  by  requiring  L?  to   correspond  to  a  third  order  formula  of  minimum 
stepnurnber: 

Lh y(t)  =  -  B{iry(t)  +  3y(t " h)  ■  i y(t  - 2h)  +  g-^*  -  3h>  +  hy,(t)}. 

The  expression  in  braces   corresponds   to  the  third  order  backward 
differentiation  formula.     This  forces 


Lj  y(t)  -  -  y(t)  +  y(t  -  h)  +  If  hy'(t)  +  |hy'(t  -  h)   -  §-y'(t  -  2h) 

+  i^y'(t  -  3h)     , 

which  corresponds   to  a  linear  combination  of  the  third  order  Adams -Bash  forth 

2         (3)       13         (3) 
and  Adams-Moulton  formulas,   namely  y=-  ABFV    '  +  jr  AMP   '.     The  pure  order 

blended  formula  is  defined  to  be 

{- y   +y    n  +  ^hf  +|hf   ,4f   o  +  i4f    o> 

-'n       Jn-^       36       n       6       n-1       4     n-2       18     n-3 

-^(4^-y    +3y     ,-|y    o  +  i-y    -,  +  hf  }  =  0. 
6     6     Jn         •'n-l       2  Jn-2       3     n-3  n 

Again   this  has  the  same  absolute  stability  region  as   the  third  order 

Enright  formula.     Note,   however,   that  this  is  a  3-step  formula  whereas  the 

mixed  order  approach  yields  only  a  2-step  formula. 

In  the  general    case  where  L,    is  a  p-th  order  m-derivative  k-step 

operator,  we  write 

i      .J  il     d,2  ,    .  m-1     d         ,  m 

Lh  "  Lh  +  h  dt  Lh  +  •  •  •  +  h       TmT  Lh 

dt 
where 

LJh  y(t)   =   -   c^  y(t)    -    ...    -   Jq  y(t   -  qh)   +  h    0fj  y'(t)    +    ... 

+  h  6^  y'(t  -  qh) 

and  q  =  max{k,   p}.     The  decomposition  is   unique   if  L   ,  L  ",...,  L     are 
chosen   (in  that  order)  so  that 

a^  =  0       ,       p  <  i   *  q   , 


and 


LJh  y(t)  =  0(hp+1)   . 


TO  m 

Having  obtained  L,  L,  ,   ...,  L,  we  define  the  associated  pure  order  blended 

formula  as  in   (2.1 ). 

The  tnright  formulas  are  attractive  because  of  their  ideal   stability 

behavior  near  hx  =  0  and  near  hx  =  °°,   and  so  it  would  be  of  interest  to 

examine  the  corresponding  blended  formulas  in  the  general   case.     Enright's 

(k  +  l)-st  order  formula 

k-1  ? 

yn   "Vl   =h   ^/li   Vi   +h     620fn 


corresponds  to  the  mixed  order  blended  formula  given  by 

* 
k 

where 


{AMF^k+1h  +  k  y*  hJ{BDF^kh 


AMF^K    u:         y-y1=hZ3-f-, 
Jn       Jn-\  .  q     l     n-i    ■ 

BDF^M      :  Z     a.  y     .   =  h  f     . 

i=0     '     n_1  n 

This   relationship  between  the  three  formulas  is  easily  verified  once  we 
know  that  6k  =   (-l)k  y£  and  ak  =   (-l)k_1/k. 

Note  that  the  order  of  the  AT4F^k+1^    -  BDF^  blend  is  not  affected 
by  changing  the  coefficient  k  y.*j   and  hence  we  consider  the  formula 

{AMF(k+1)}   -  y^   hJ{BDF(k)} 

(k)  (k) 

where  y         is  a  free  parameter.      In  particular  y         could  be  chosen  to 

maximize  the  angle  a  for  which  the  formula  is  A(a)-stable.      (A  formula 

is  A(a)-s table  if  e^jery  numerical    solution  of  y'    =  xy  tends  to  zero  for 

|Arg(-x)|    <  a.)     Values   for  y         numerically  determined  in  this  way  appear 


(k)         (k) 
in  Table  I  in.  the  column  yv    '  =  YqD1.     The  largest  angle  a  for  which  ea 


ch 


10 

formula  is  A(a)-stable  is  also  given.     A(a)-stable  formulas  up  to  order  12 
were  obtained.     The  second,   third,   and  fourth  order  formulas  are  A-stable 
for  a  range  of  values  of  y       »  y       »  and  y       .     The  actual    choice  of  these 
three  parameters  is   delayed  until   §4. 

(k) 
Table   I.     Values   for  yv    '   and  Corresponding  a-angles 


a  for 

a  for 

Stepnumber 
k 

Order 
k  +  1 

- k  Yk 

(k) 

Yopt 

(k)           l     * 
yv  '  =  -  k  Yk 

(k)  ..     (k) 

Y         "  Yopt 

1 

2 

.5 

[  0,+-) 

90° 

90° 

2 

3 

.1666667 

[.1250000,+°°) 

90° 

90° 

3 

4 

.125 

[.1218908, 
.6837917] 

90° 

90° 

4 

5 

.1055556 

.1284997 

86.5° 

89.4° 

5 

6 

.09375 

.1087264 

79.6° 

87.0° 

6 

7 

.08561508 

.09625961 

72.5° 

82.9° 

7 

8 

.07957176 

.08754865 

61.5° 

77.4° 

8 

9 

.07485229 

.08105623 

40.2° 

70.2° 

9 

10 

.07103299 

.07599874 

not  A(0)-s table 

60.7° 

10 

11 

.06785850 

.07192936 

not  A(0)-stable 

47.6° 

11 

12 

.06516462 

.06857227 

not  A(0)-stable 

28.7° 

The  second  last  column  was  obtained  from  measurements  of  Figure  3  in 
[4,   p.    21]. 


The  angle  a  is  only  one  of  a  number  of  parameters  which  have  been 
proposed  for  measuring  the  extent  of  the  stability  region  R.     But  it  is 


11 

probably  the  best  such  measure,  especially  for  methods  with  automatic 
stepsize  selection.     When  such  methods  are  applied  to  the  equation 
y'   =  Ay  +  g(t)  with  A  complex,  the  integration  starts  out  with  hA  near  the 
origin,  and  as  the  integration  proceeds,  the  stepsize  h  increases  and  the 
value  hA  moves  away  from  the  origin.     However  if  hA  approaches  the  boundary 
of  R,   this  would  be  detected  by  the  error  estimator,  and  any  further  move- 
ment of  hA  would  be  prevented.     The  implication  of  this  is   that  not  all   of 
R  is   "used",  but  only  that  portion  which  can  be  reached  from  the  origin  by 
rays  lying  entirely  inside  R.     The  useful   part  of  R  can  be  described 
mathematically  as  the  set 

{y  e  R:      t  p  e  R  for  0  <  t  ^   1}, 
and  for  stiffly  stable  formulas  this  region  is  much  like  a  wedge,  which  is 
best  described  by  the  angle  a.     This  parameter  also  serves  as  a  good 
indicator  of  the  problem  class  for  which  the  formula  is  suitable,  namely 
those  problems  for  which  the  "large"  eigenvalues  A  of  the  Jacobian  satisfy 
| Arg(-x) |    <  a.     Ideally  we  would  like  to  have  o  =  5-  so  that  the  formula 
would  be  suitable  for  all   realistic  problems. 

It  may  be  helpful   to  present  the  original  motivation  for  the 

(k+1)  (k) 

AMFV  -  BDFV       blend.     We  begin  by  considering  a  scalar  equation 

y1    =   f(t,  y).      If  the  stepsize  h  were  determined  only  by  accuracy 

considerations,   then  the  size  of  -h  f    would  be  a  good  measure  of  stiffness. 

If  -h  f    were  small,  then  the  AMP  + ''  would  be  a  good  k-step  formula  to 

use,  and  if  -h  f    were  large  then  the  BDFV    ;  would  be  a  good  choice.     And 

so  it  would  seem  that  a  weighted  sum 

{AMF(k+1)}  +  y(k)(-h   f  ){BDF(k)} 


12 
would  give  us  the  best  of  both   formulas.     For  small    -h  f     this  blend  is 

y 

nearly  the  AMPk+1\   and  for  large  -h   f     it  is  nearly  the  BDPk^. 

We  now  extend  this  idea  to  a  system.      Consider  a  change  of 
coordinates  which  decouples  the  system  at  (t,  y)  =  (t  ,  y  ).     Such  a 
coordinate   transformation  can  be  accomplished  by  a  matrix  T  which 
diagonalizes  f  (t  ,  y  );   that  is,    f"     fv(t  ,  Yn)T=  A  where  a  is  diagonal. 
This   is   almost  always  possible.     Letting  z  =  T~     y,   the  system  becomes 

z'  =  g(t,  z)  =  T"1   f(t,  Tz). 
Since  this  system  is  locally  decoupled,  we  can  use  -h  a11  as  a  measure  of 
the  stiffness  of  the   i-th  differential   equation.     Each  equation  of  the 
system  can  be  integrated  by  a  formula  which  is  suitable  for  that  equation: 

Transforming  this  back  to  the  original   coordinate  system  gives 

K  +  Vl  +  h  J0  6J   W  +  T{k)(-h  V("j0  'i  yn-.i  +  h  fn>  "  °- 

which  is   the  blended  formula. 

The   remainder  of  the  paper  is   restricted  to  a  discussion  of  the 

AMF^k+1^   -  BDF^   blend  with  Y^   =  y^i- 

'opt 
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3.     THE  NORDSIECK  REPRESENTATION 

In  the  Nordsieck  representation  we  work  with  an  approximation 
^  =  [yn>  h  y^,   ...,  hk  y^/k!]T  to  the  scaled  derivatives  [y(tn), 
h  y'(t  ),   . ..,  hk  y^k^(tn)/k!]T  of  the  solution.     A  linear  Nordsieck  formula 
can  be  thought  of  as  consisting  of  a  prediction 

^  =  P  Vl 

followed  by  a  correction 

*n  =  En  +  *■  An 

where  a  is  chosen  so  that  y1  =  f(t  ,  v.,);  that  is,  a  is  defined  implicitly 
n  •'n     n   n  n  v  J 


by 


An  =  ju^h  f(t  ,  pn  +  £n  a J  -  h  p'} 


Here  p^  =  [pn>  h  p^,  . . . ,  hk  p^/k!]1,  £  =  [£Q,  a} £k]T,  and  P  is  the 

Pascal  triangle  matrix  given  by 

Pid  =  (  ]  )       ,   0  «  i,  j  $   k. 

Nordsieck  formulas  are  closely  related  to  multistep  formulas.  For 

every   linear  k-step  formula  of  order  >,   k  there  is  a  unique  (k  +  l)-value 

linear  Nordsieck  formula  such  that  the  values  y  ,  y     -,,...,  y  ■  generated 

by  the  Nordsieck  formula  from  v  .  satisfy  the  difference  equation  of  the 

multistep  formula.  Let  the  multistep  formula  be 

p(E)  y  ,  =  c(E)  h  f  . 
x  .  Jn-k        n-k 

where  E  is  the  shift  operator, 

p(0  =  Z     a.  £   ,   and  aU)   =     E  B  5k_1  . 
1=0  ]  i=0  ^ 

The  corrector  vector  _£  of  the  "equivalent"  Nordsieck  formula  satisfies 
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p(c)  =  u   -  Dk+1  &f(5  I  -  P)'1  L  (3-D 

and 

oU)   =  (E  -  Dk+1  ej(e  I  -  P)"1  *  (3.2) 

where  e,  =  [0,  1 ,  0,  . . . ,  0]  and  e«  =  [1,  0,  0,  ...,0].  It  can  be  shown 
from  these  equations  that 

l-   =   aQ   and  zQ   =  60  .  (3.3) 

Proofs  of  these  relationships  will  appear  in  [19]. 

(k)     (k) 
Let  a_v  '   and  tr  '  denote  the  corrector  vector  corresponding  to  the 

k-step  AMPk+1)  and  BDPk'  respectively.  It  can  be  shown  [19]  that 

a(k)  =  .k;]  (.dj.iw) 


and 


b(k)  u   .lr  k+1  1        •  m    Q(1)k 

Dj     k!  L  j+1  J   '   J   UU'K' 


where  the  brackets  denote  Stirling  numbers  of  the  first  kind  (see  e.g.  Knuth 
[13,  p.  66]).  Values  of  a^^  and  b/k^  are  given  in  Appendix  A  for  k  =  1(1)11. 

Consider  the  Nordsieck  formula  with  £  =  a  •  y  h  J  b ,  where  we  have 
suppressed  the  superscripts  (k).  Applying  (3.1)  and  (3.2)  gives 
p(0  =  pAMU)  -  y  h  J  pBDU)  and  aU)   =  aAM(?)  -  y   h  J  aBD(s).  Therefore 
the  values  computed  by  this  Nordsieck  formula  satisfy 

P   (E)  yn_k  -ha   (E)  fp_k  -  y  h  J{p   (E)  yp_k  -ha   (E)  fp_k>  =  0  , 

(k+1)      (k) 
which  is  the  difference  equation  for  an  AMFV   '  -  BDFX  '   blend. 
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In  an  actual   implementation  J  might  be  reevaluated  every  few  steps; 
in  other  words, 

where  i     =  a  -  y  h  J     b.     And  a  blended  Nordsieck  formula  with  a  changing 
-n      —      '         n  —  3     3 

J  is  not  equivalent  to  a  blended  multistep  formula  with  a  changing  J.     For 
example,  the  blended  3-value  Nordsieck  formula  generates  values  y  ,  y, , 
y  _p  which  satisfy 

<!  +  3y  h  VX  -  Vl   -  fl  fn  -  tK-1 


(3.5) 


-Yhyfyn-fVl  -hfn-|W] 
♦  0  ♦  3r  h -W^C  jfc  V,  +  t|  V2 

This  difference  equation  can  be  derived  as  follows.  Premultiply  (3.4)  by 
e_,  to  get 

M;  =  ejpVl  +  ilnAn  , 

and  substitute  this  back  into  (3.4): 

where  I     =  aZ     i.     From  this  equation  it  follows  that 
— n         In  — n  n 

yn.,  =  4(I  -  Vi «?'  p  Vz +  4  !„_!  h  y;.i 

and 

Vn  ■  ej(l  -  i„  e{)  P(l-Vl4>  PV2 

+  e^(I  -!„*[>  P  Vl   hj,n-l  +  §0^hyn   ' 
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which  is  a  system  of  two  equations  relating  y.  y  , ,  y  „  hy',  h  y1  . 
h  y;_2,  and  h2  y^/2. 
yields  equation  (3.5). 


2  2 

h  y'  os  and  h  y"  0/2.  Eliminating  h  y"  0/2  and  substituting  f  for  y1 
n-d        n-£  n-£  n     n 
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4.  THE  CORRECTOR  ITERATION 

The  implicit  equation  satisfied  by  a  is 

h  p1  +  a,  A  -  h  f(t,  p  +  &0  a)  =  0 

where  for  convenience  we  suppress  dependence  on  n.  The  Newton  iteration  for 
this  equation  is 

A(0)  =  °» 

A(m+1)   =  A(m)  +  ^1   -  h(fy)(m)  ^0]_1{h  f(m)  "   (h  P'   +  *1  A(m))} 
where  for  a  function  the  subscript  (m)  denotes  evaluation  at  (t,  p  +  £Q  A/    v). 
This  iteration  is  the  same  for  both  the  multistep  form  and  the  Nordsieck  form 
of  the  formula.     Convergence  of  this  iteration  is  seldom  a  problem,  and  so  a 
considerable  saving  in  time  results  if  a  recent  Jacobian  approximation  J  is 
used  in  place  of  (fy)(m): 

A(m+1)  =  A(m)  *CVhJ  V"1{h  f(m)   "    (h  P'   +  £1  A(m))}'     <4'« 
From  i  =  a  -y  h  J  b  and  from  (3.3)  we  have  that 

£,   -  h  J  iQ  =  (1   -  y  a  h  J)  -  h  J(b  -  y  h  J) 

BD         AM 
where  a  =  an  and  3  =  $Q  .  This  means  that  (4.1)  involves  the  LU  decomposi- 
tion of  a  quadratic  polynomial  in  h  J.  But  forming  matrix  products  can  be 
expensive,  and  this  is  particularly  true  if  sparse  matrix  techniques  are  being 
used  because  matrix  multiplications  reduce  sparseness. 

One  way  of  avoiding  matrix  multiplications  is  suggested  by 
Willoughby  [6  ,  p.  102],  and  it  is  applicable  as  long  as  the  roots  of 
1  -  (y  a  +  e)z  +  y  z  are  complex  conjugates.  This  technique  calls  for  an 
LU  decomposition  in  complex  arithmetic,  which  has  the  disadvantage  of  being 
four  times  slower  and  requiring  twice  as  much  storage  compared  to  real 
arithmetic. 
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A  second  possibility  is  to  choose  the  parameter  y  so  that 

2  2 

1  -  (y  a  +  3)2  +  y  z  is  a  perfect  square  (1  -  c  z)  .  Then  multiplication 

-2 

by  (1  -  c  h  J)   can  be  accomplished  by  decomposing  1  -  c  h  J  and  performing 

two  backsolves.  Enright  [6]  derives  a  second  set  of  formulas  based  on 
this  idea.  This  approach,  however,  reduces  the  number  of  free  parameters, 
and  so  one  order  of  accuracy  must  be  sacrificed. 

There  is  a  third  way  of  avoiding  matrix  multiplications,  and  that 

2  2 

is  to  approximate  1  -  (y  a  +  b)z  +  y  z  by  a  perfect  square  (1  -  c  z)  so 

that  instead  of  (4.1)  we  use 

<Vl)  "  %)  +  V   -   c  h  J]"2{h  f{||)  -  (h  p-  ♦  *,  A(m))}.   (4.2) 

This  is  the  iteration  which  is  used  in  our  implementation  of  the  blended 
formulas. 

Recall  that  for  k  =  1,  2,  3  there  is  some  freedom  in  the  choice  of 

2 
y.  We  use  this  freedom  to  make  1  -  (y  a  +  b)z  +  y  z  as  close  as  possible 

to  a  perfect  square.  There  are  two  choices  of  y  which  yield  a  perfect 

square: 

y(1)=f±/2  , 

y(2)  =  Jl+  2  £ 

18  *  9  /D  ' 

(3)  =  189  +  J8  ^ 
Y     484  -  121  °  ' 

The  smaller  value  of  y  in  each  case  results  in  a  smaller  truncation  error, 

and  so  we  choose  the  y's  as  close  as  possible  to  the  smaller  value  subject 

to  the  constraints  imposed  by  A-stability  (see  Table  I): 

y(1)  =  |  -  /2   =  .08578644  , 

Y^  =  .1250000  , 

y(3^  =  .1218908  . 
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The  iteration  (4.2)  is  significantly  different  from  (4.1)  in  that 
it  converges  linearly  for  the  system  y'  =  J  y  whereas  (4.1)  converges  after 
one  iteration  for  this  system.  Thus  it  is  important  to  choose  c  so  that 
convergence  is  rapid  for  this  special  case.  In  practice  the  eigenvalues 
of  h  J  have  either  a  small  modulus  or  a  negative  real  part.  However,  it  is 
only  those  eigenvalues  with  negative  real  part  which  might  cause  convergence 
difficulties.  Hence  the  value  of  c  is  chosen  by  considering  the  test 
equation  y'  =  \  y,  Re  A  $  0.  For  this  problem  the  iteration  given  by  (4.2) 
becomes 

A(m+1)   =  A(m)  +H   -  ch  A]"2{h  X[p+  (6  -yh  A)A(m)] 

-    [h    p'    +    (1    -   Y   h    A   a)A(m)]}    . 

Letting  y  =  h  A  and  simplifying  gives 

A(m+1)  =  [1  "  c  ^~2{^  a  +  e  "  2c)y  +  (°2  "  y)y2}A(m) 

+  [1  "  C  y]"2  h(X  P  -  p')  . 

The  rate  of  converaence  depends  on  the  coefficient  of  Af    \,  and  so  we  choose 

(m) 

c  to  minimize 

♦  (c)   ■  Rem^  0    |[1   -  c  y]_2{(Y  a  +  3   -  2c)y  +  (c2  -  y)u2>|    . 

Clearly  we  must  choose  c  >  0,  and  so  the  expression  inside  the  vertical 
bars  is  analytic  for  Re  y  $  0.     By  the  maximum  modulus  theorem  the  maximum 
is  attained  for  Re  y  =  0,  and  so 

max  |(t  «+B   -  2c)1   M-   (c2  -T)M2| 

-°°<a)<+°°  |1    -   c   l    bl| 

Equivalently 
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Mc)f  ,  max  »'[(T  .  >  b  -  2c)2+  (c2  -tfJl 

\_\    +  C     w   J 
Differentiation  with  respect  to  u  yields  extrema  at 

J  =  0, 
2 


iii 


and 


2 

OJ     = 


(y  a  +  3   -  2c) 


c2(y  a  +  6   -   2c)2   -   2(c2   -  y) 


2         if  2(c2  -  Y)2  <  c2(y  a  +  6   -  2c)2; 


therefore 


V  -  ^ 


if         2(c2  -  Y)2  *  c2(Y  a  +  6  -  2c)2   , 


Mc)f  -  ( 


max  { 


(c2  -  v)2 


v« 


(y  a  +  3   -  2c)' 


4[c2(y  a  +  6    -   2c)2   -    (c2   -  y)2] 


}         otherwise 


The  value  of  c  which  minimizes  41(c)  was  determined  numerically  for 
k  =  1(1)11.     These  values  appear  in  Table  II. 
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Table  II.  Values  for  c 


(k) 


Stepnumber 
k 

(k) 

*(c<k)) 

1 

.2928932 

0 

2 

.3374973 

.1184662 

3 

.3335427 

.1161818 

4 

.3427329 

.1139836 

5 

.3169058 

.0995555 

6 

.2992971 

.0894459 

7 

.2862392 

.0819169 

8 

.2760327 

.0760633 

9 

.2677630 

.0713660 

10 

.2608834 

.0675036 

11 

.2550426 

.0642651 

From  Table  II  we  see  that  even  in  the  worst  case,  x  pure  imaginary 
and  k  =  2,  the  asymptotic  error  coefficient  is  only  .118,  which  means  that 
each  iteration  is  good  for  an  additional  0.93  digits. 
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5.  ERROR  ESTIMATION  AND  CONTROL 

The  purpose  of  estimating  the  local  error  is  to  assist  in  the 
selection  of  stepsize  and  order  in  such  a  way  that 

(i)   the  least  possible  work  is  done  for  the  accuracy  attained. 

(ii)  the  global  error  behaves  in  a  reasonable  manner  as  a  function 
of  the  local  error  tolerance. 
Hence  the  estimation  of  the  local  error  is  not  in  itself  of  much  importance 

The  theory  of  error  estimation  is  not  well  developed.  Usually  an 

attempt  is  made  to  devise  a  local  error  estimate  which  is  asymptotically 

equal  to  the  local  error  under  the  assumption  that  the  order  and  stepsize 

do  not  vary.  This  approach  must  be  used  cautiously.  For  example  the 

(strongly  A-stable)  formula  y     =  y     ,  +   .501  h  f  .,  +  .499  h  f  has  the 
v    3  J  '  Jr\      -'n-l         n+1         n 

asymptotically  correct  error  estimate  .001  h(f  .,  -  f  ).  However,  this 
will  often  give  a  \iery   poor  estimate,  simply  because  the  leading  term  in 
the  asymptotic  expansion  of  the  local  error  is  not  the  dominant  source  of 
error  unless  h  is  \ievy   small.  A  more  realistic  approach  would  be  to  make 
use  of  a  good  bound   on  the  truncation  error.  This  approach  would  yield 

the  estimate  .250001  hi  If  -  f  ,1  for  our  example.  It  is  worth  noting 

1  '  n    n- 1 '  '  r 

that  the  highly  rated  code  DE  of  Shampine  and  Gordon  [17]  uses  error 

estimates  that  are  not  asymptotically  correct;  for  example,  the  estimate 

.5  h(f  -  f  , )  is  used  for  the  trapezoid  formula. 
v  n    n-1  K 

Our  implementation  of  the  blended  formulas  uses  an  error  bound 
as  the  basis  for  the  error  estimate.  In  order  to  derive  the  error  bound, 
we  first  split  the  AMF^k+1^  into  a  sum  of  the  AMF^  and  the  ABF^: 

AMF<k+1>  =^AMF<k>  +^ULAABF(k) 
Yk-1  Yk-1 
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where  y.    is  given  in  Henrici  [11,  p.  193]  and  Gear  [10,  p.  108].  This 
equality  follows  from  the  fact  that  eQ  =  y,  for  the  AMP   '.  From  Henrici 
[11,  equation  (5-22)]  we  have  that  y.    -  y*-.   =  y* ,  and  hence 

Yk-i      Yk-i 

If  we  consider  the  AMF-BDF  blend  applied  to  a  scalar  equation,  we  get  that 
the  (normalized)  truncation  error  is 

hw(1  .  Y  h  x)-igJiy<M>(t)  .  IkJly(W)(T.,  .  rU,<w){t.,]  . 

Yk-1  Yk-1  K   ' 

And  since 


Yl  YP 


we  have  the  bound 


2—-^  *  k-TT   for    k  =  1(1)]1 
Yk-1     K   ' 


hk+1      max     |v(k+l)m 
k  +  1  t    <  t  <  t   "     *L' 


on  the  truncation  error  when  \   is  real  and  negative.  Hence  we  use 

k  +  1  I  ,Jfn    yn-l  I  ' 
as  our  error  estimate. 
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6.      COMPARISON  WITH  THE  BACKWARD  DIFFERENTIATION  FORMULAS 

We  begin  by  examining  the  changes  that  must  be  made  to  a  BDF 
algorithm  in  order  to  convert  it  into  an  algorithm  for  the  AMF-BDF  blend 
The  prediction  step  is   identical   for  the  two  formulas,  and  it  is  only  in 
the  correction  step  that  there  are  differences.     The  BDF's  have  a 
correction  step  something  like 

A  «-  0 

for  m  «-  1 ,  2,  3 
rf  +  f(t,  y) 


if  desirable  then 


fw  «-  b,  -  h  x  f  (t,  y) 

v 


do  / 


decompose  W 
solve  W  x  d  =  h  x  f  -  hy'   for  d 

y  +■  y  +  d 

hy'  «-  hy'   +  b]  *  d 

A   •*-   A   +   d 

rf  d  small   enough   then  £p_  to  testerror 
deal  with  iteration  convergence  failure 
testerror:     rf  A/(k+l)   too  large  then  go  to  tryagain 


v. 


for  j  *-  2(1  )k  do  hjyuy/j!  *■  hJy^JVJ!  +  b.    x  a 


,Jy(J) 

whereas  for  the  AMF-BDF  blend  we  have 
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r 


A  +■  0;  A  ■«-  0 
f  o  r  m  «-  1 ,  2 ,  3 
>  *  f(t,  y) 

'""cH  «-  c  x  h 
if  desirable  then  /  W  <-  1    -   cH  x  f  (t,  y) 

^decompose  W 
solve  W  x  d  =  h  x  f  -  hy'   for  d 
do_  /  solve  W  x  d  =  d  for  d 
d^Yxh><   (d-  d)/cH 
y  *-  y  +  an  x  d  +  d 
hy'  ^  hy'   +  d  +  b1   x  d 
A^A  +  d;A«-A  +  d 

rf  d  +  d  small  enough  then  go  to  testerror 
deal  with  iteration  convergence  failure 
testerror:     rf   (a  +  A)/(k  +1)   too  large  then  go  to  tryagain 

for  J  +  2(1  )k  do  hVJ"Vj!  «■  hV     AJ1  +  a.   *  A  +  b.  x  a 
In  both  cases  one  function  evaluation  is  required  for  each   iteration. 
However  the  corrector  overhead  for  the  AMF-BDF  blend  is  twice  as  great  as  for 
the  BDF.     The  overhead  is  defined  to  be  any  computation  other  than  function 
evaluations. 

A  significant  portion  of  the  overhead  of  a  stiff  system  solver 
often  consists  of  LU  decompositions.     We  would  not  expect  there  to  be  a 
difference  between  the  two  formulas  in  the  number  of  decompositions  if  we 
were  to  use  the  same  type  of  iteration  in  both  cases.     For  the  BDF  we  use 

A,„.,x    =    A, _,    +    [£,    -    H    J    £n]_1{...} 


V." 


(m+1)   ~  "(m) 


1 


'0- 
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where  H  is  the  size  of  the  step  when  J  was  last  evaluated.  The  same  type  of 
iteration  for  the  blended  formula  would  be 

A(m+1)  =  A(m)  +  [(1  "  y  «  h  J)  -  H  J(g  -  Y  h  J)]"1!...} 
However  what  we  actually  use  is 

Vl)  =A(m)  +  [1  -cH  J]-2(...}. 
This  second  iteration  might  be  expected  to  fail  more  often  than  the  first, 
but  it  is  difficult  to  analytically  compare  the  two  iterations  even  for  the 
equation  y'  =  Ay,  and  therefore  we  have  to  rely  mainly  on  empirical  evidence. 

If  the  preceding  discussion  is  generalized  to  blends  of  m  formulas,  it 
can  be  seen  that  the  correction  iteration  overhead  is  proportional  to  m. 
Hence  blends  of  m  formulas  are  competitive  only  if  they  require  significantly 
fewer  function  evaluations  to  achieve  a  given  accuracy.  How  many  fewer 
depends  on  the  ratio  of  the  cost  of  a  function  evaluation  to  the  cost  of  a 
backsolve.  For  problems  with  very   simple  right  hand  sides  f(t,  y) ,  the  cost 
of  a  function  evaluation  might  be  roughly  equal  to  the  cost  of  a  backsolve 
and  in  this  case  it  would  not  pay  to  use  a  blend  of  m  formulas  unless 
the  number  of  function  evaluations  compared  to  that  of  the  BDF's  is  smaller 
by  a  factor  of  2/(m  +1).  Here  we  are  assuming  that  the  number  of 
decompositions  is  independent  of  m,  which  is  fairly  reasonable  since  it  is 
primarily  the  heuristics  used  by  the  program  that  determine  the  number  of 
decompositions. 

We  have  compared  the  efficiency  of  the  blended  formulas  to  that  of 
the  BDF's  on  several  test  problems.  In  order  to  make  the  comparison  fair, 
the  blended  formulas  were  implanted  into  a  well-known  BDF  code  (Gear's 
DIFSUB  [9])  in  such  a  way  that  only  the  formula-dependent  part  of  the  code 
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was  changed.  And  therefore  the  code  was  not  tuned  for  the  blended  formulas. 
There  was  one  modification  made  to  the  codes  and  that  was  to  replace  the 
call  to  MATINV  by  calls  to  DECOMP  and  SOLVE  [8].  A  listing  of  the  source 

ft 

code  appears  in  Appendix  B. 

In  order  to  facilitate  testing,  only  problems  with  known  analytic 

solutions  were  used.  The  stepsize  h  was  initially  set  to  the  local  error 

tolerance  e,  and  on  subsequent  steps  it  was  set  to  min{h,  tf  -  t}  where  h 

is  the  stepsize  recommended  by  the  method  and  [0,  tf]  is  the  interval  of 

integration.  The  parameter  h    was  set  to  tf  -  t  and  h  .  was  (mistakenly) 
3  max  t        mm 

1  -13 

set  to  j  u  tf  where  u  is  the  unit  roundoff  error,  16   ,  of  the  IBM  360. 

The  weight  vector  w  was  initialized  to  all  ones.  The  Jacobian  f  (t,  y)  was 

evaluated  numerically..  The  relative  (global)  error  was  defined  by 

i    i  /  \ 

max  /  J  ,yn  "  y  ^n\2xl/2 

-  '  OsnsN  {.   ,  *     i     ;  ; 
i=l 

where  w  =  max{l,  yn ,  y, ,  ...,  y  }. 

Tables  III  through  VII  give  the  results  of  comparing  DIFSUB  with 

blended  DIFSUB  for  five  different  types  of  problems.  Each  problem  was 

-2    -3        -10 
integrated  for  nine  different  error  tolerances,  e  =  10  ,10  ,  ...,  10 

For  each  integration  the  following  statistics  were  collected: 

max  order      -  the  order  of  the  highest  order  formula  selected  by  the 

code. 

steps         -  the  total  number  of  steps  taken. 

func  evals     -  the  total  number  of  function  calls  (including  those 

needed  to  approximate  the  Jacobian). 

back  solves    -  the  total  number  of  calls  to  the  subroutine  SOLVE, 

which  solves  a  linear  system  that  is  in  factored  form. 

LU  decomps     -  the  total  number  of  calls  to  the  subroutine  DECOMP, 
which  performs  an  LU  factorization  on  a  matrix. 
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accurate  digits  -  the  negative  of  the  base  10  logarithm  of  the  relative 

(global)  error,  which  was  defined  previously. 

time  -  machine  time  in  centiseconds  on  the  IBM  360  model  75 

at  the  University  of  Illinois  at  Urbana-Champaign. 
These  timings  are  not  completely  reliable  due  to 
multiprogramming. 

It  should  be  stressed  that  the  intent  is  to  compare  two  classes  of 

formulas   and  not  two  methods  or  two  programs.  Also,  only  one  aspect  of 

the  relative  utility  of  the  two  formulas  is  being  compared,  namely  machine 

time  versus  global  accuracy. 

Problem  I.     The  following  is  a  linear  problem  whose  Jacobian  matrix  has 
eigenvalues  -.1,  -50,  -120: 


-.1 

-49.9 

0 

0 

-50 

0 

0 

70 

-120 

y(o)  ■ 


This  was  integrated  on  [0,  15] 
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Table  III.  Numerical  Results  for  Problem  1 


max 

steps 

func 

back 

LU 

accurate 

time 

£ 

order 

evals 

solves 

decomps 

digits 

DIFSUB 

ID"2 

3 

29 

104 

76 

9 

1.9 

30 

_3 

10 

5 

49 

145 

108 

12 

2.9 

38 

io-4 

4 

70 

202 

171 

10 

3.9 

51 

10"5 

6 

106 

286 

240 

15 

5.0 

75 

10"6 

5 

193 

508 

453 

18 

5.7 

128 

IO"7 

6 

176 

474 

422 

17 

6.7 

128 

10"8 

6 

204 

551 

505 

15 

7.5 

143 

IO"9 

6 

270 

771 

719 

17 

8.4 

185 

IO"10 

6 

353 

1024 

969 

18 

9.3 

251 

bier 

ided  DIFSUB 

IO"2 

4 

30 

101 

146 

9 

3.2 

33 

10   J 

5 

43 

141 

220 

10 

3.9 

45 

IO"4 

6 

69 

195 

316 

12 

4.6 

71 

IO"5 

7 

89 

238 

396 

13 

5.4 

86 

IO"6 

8 

120 

308 

524 

15 

6.5 

no 

IO"7 

11 

139 

410 

662 

26 

7.4 

155 

IO"8 

9 

169 

443 

782 

17 

8.5 

165 

IO"9 

11 

209 

540 

952 

21 

9.3 

228 

IO"10 

11 

234 

595 

1056 

22 

10.4 

261 

An  examination  of  the  last  two  columns  reveals  that  there  is  a 
slight  though  definite  improvement  in  performance  when  the  blended 
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formulas  are  substituted  for  the  BDF's.  The  reason  that  blended  DIFSUB 
does  more  LU  decompositions  might  be  due  to  the  fact  that  fresh  decomposi- 
tions are  performed  whenever  the  order  is  changed  but  not  when  the  stepsize 
is  changed. 

Problem  II,     The  following  is  problem  12  in  Krogh's  [14]  set  of  test 
problems;  it  is  also  used  by  Gear  [10,  p.  218]: 

-  800z^  +  {z? 


y'  =  u 


-lOOOz1  +  (z1' 


.2  .  ,_2 


where  z  =  Uy  and 


10z3  +  (z3 

4   /  4 
-.001zH  +  (z* 


-1 


U=I 


1 


y(o)  = 


-l 
-l 
-l 


l 
-l 


l   -i 


This  was  integrated  over  [0,1000].  The  eigenvalues  of  the  Jacobian  matrix 
are  equal  to  -1002,  -802,  8,  and  -2.001  at  t  =  0,  and  approach  -1000,  -800, 
-10,  and  -.001  respectively  as  t  ■*■  ». 


31 
Table  IV.  Numerical  Results  for  Problem  2 


e 

max 
order 

steps 

func 
evals 

back 
solves 

LU 
decomps 

accurate 
digits 

time 

DIFSUB 

10"2 

3 

64 

277 

184 

23 

1.7 

65 

10"3 

4 

105 

374 

273 

25 

2.7 

91 

ID'4 

5 

167 

496 

387 

27 

3.6 

131 

lO"5 

6 

250 

733 

608 

31 

4.1 

206 

10"6 

6 

284 

778 

677 

25 

5.4 

229 

lO"7 

6 

380 

1076 

915 

40 

6.4 

318 

10"8 

6 

467 

1311 

1178 

33 

7.0 

416 

lO"9 

6 

590 

1641 

1516 

31 

8.2 

586 

lO"10 

6 

768 

2083 

1954 

32 

9.1 

699 

bier 

ided  DIFSUB 

in"2 

4 

61 

276 

358 

24 

2.7 

75 

lO"3 

6 

98 

333 

488 

22 

3.7 

105 

lO"4 

7 

146 

449 

688 

26 

4.6 

160 

10  D 

7 

200 

594 

970 

27 

5.4 

231 

10"6 

9 

262 

750 

1242 

32 

6.4 

314 

lO"7 

9 

325 

900 

1502 

37 

7.3 

419 

10-8 

9 

381 

1062 

1818 

38 

8.4 

521 

lO"9 

11 

474 

1303 

2244 

45 

9.4 

734 

10"10 

11 

558 

1548 

2678 

52 

10.4 

852 

Once  again  there  is  a  slight  though  definite  improvement  when  the  blended 
formulas  are  substituted  for  the  BDF's. 
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y.  = 


y(o)  = 


Problem  ill.     This  is  problem  B5  in  the  set  of  stiff  problems  considered  by 
Enright  et  al  [7]: 

0 

0 

0 

0 

0 

-.1 

This  was  integrated  over  the  interval  [0,  20].  Because  the  Jacobian  matrix 
has  eigenvalues  -10  ±  ilOO,  -4,  -1,  -.5,  -.1,  this  problem  can  be  troublesome 
for  the  BDF's. 


-   10 

100 

0 

0 

0 

-100 

-10 

0 

0 

0 

0 

0 

-4 

0 

0 

0 

0 

0 

-1 

0 

0 

0 

0 

0 

-. 

0 

0 

0 

0 

0 
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Table  V.     Numerical   Results  for  Problem  3 


£ 

max 
order 

steps 

func 
evals 

back 
solves 

LU 
decomps 

accurate 
digits 

time 

DIFSUB 

io-2 

(4) 

(1001) 

(3002) 

(2959) 

(7) 

(1.1) 

"too 

much  work' 

'  ;  integration  s 

topped  at  t 

=  8.3 

blende 

d  DIFSUB 

ID"2 

6 

125 

493 

756 

19 

2.9 

215 

10"3 

6 

194 

691 

1152 

19 

3.7 

348 

IO"4 

7 

277 

922 

1590 

21 

4.5 

517 

10  3 

8 

368 

1196 

2114 

23 

5.4 

729 

IO'6 

9 

468 

1494 

2686 

25 

6.4 

1040 

IO"7 

11 

592 

1831 

3288 

31 

7.3 

1546 

IO"8 

11 

721 

2178 

3958 

33 

8.5 

1772 

IO"9 

12 

895 

2644 

4878 

34 

9.4 

2369 

IO"10 

(12) 

(1001) 

(2917) 

(5580) 

(21) 

(10.3) 

"too 

much  work' 

'  ;  integ 

ration  s 

topped  at  t 

=  5.1 

There  is  a  remarkable  difference  between  the  abilities  of  the  two  types  of 
formulas  to  solve  this  problem. 

Problem  IV.     The  following  is  problem  13  in  Krogh's  set  of  problems;  it  is 
also  used  by  Enright  et  al   [7,  problem  E4]: 
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lOz1  +  10z2  +  \  (z1)2  -  \  (z2)2 

1     2    12 
-lOz1  +  lOz  +  z1  z 

-lOOOz3  +  (z3)2 

L   -01z4+(z4)2 


where  z  =  Uy  and 


-1 


1 

1 

1 

-1 

1 

1 

1 

-1 

1 

1 

1 

-1 

y(o)  = 


o 

-2 

-1 
-1 


u  =  i 


1 

This  was  integrated  over  [0,1000].  The  eigenvalues  of  the  Jacobian  matrix 
are  equal  to  8  ±  lOi,  -1002,  and  -2.01  at  t  =  0  and  approach  -10  ±  lOi, 
-100,  and  -.01  as  t  ->  ». 
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Table  VI.     Numerical   Results  for  Problem  4 


E 

max 
order 

steps 

func 
evals 

back 
solves 

LU 
decomps 

accurate 
digits 

time 

DIFSIJB 

io-2 

3 

79 

388 

263 

31 

1.3 

78 

10-3 

4 

137 

555 

402 

38 

2.0 

116 

IO"4 

5 

289 

836 

715 

30 

3.3 

249 

IO"5 

6 

268 

802 

677 

31 

4.2 

229 

IO"6 

(6) 

(1001) 

(2777) 

(2652) 

(31) 

(5.0) 

"too 

much  work' 

"  ;  integration  s 

topped  at  t 

=  150 

blende 

d  DIFSUB 

IO"2 

4 

79 

347 

484 

26 

2.6 

95 

IO"3 

5 

127 

463 

708 

27 

3.5 

140 

IO"4 

7 

170 

558 

890 

28 

4.7 

191 

-5 
10  ° 

7 

244 

750 

1258 

30 

5.4 

266 

IO"6 

7 

325 

952 

1646 

32 

6.4 

376 

IO"7 

9 

399 

1165 

2008 

40 

7.4 

536 

IO"8 

11 

483 

1367 

2364 

46 

8.3 

724 

IO"9 

(11) 

(1001) 

(5511) 

(6708) 

(539) 

(9.2) 

1389 

"too 

much  work1 

1  ;  integ 

ration  s 

topped  at  t 

=  977 

For  this  problem  the  blended  formulas  performed  significantly  better. 

Problem  v.     The  following  problem  (problem  8  of  Krogh's  set  of  test  problems) 
tests   the  integrator's  ability  to  solve  nonstiff  problems: 
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y    ■ 


1,,   K2       ,  2v2v-3/2 

-y  ((y  )    +  (y  )  ) 

2,,   K2       /   2v2N-3/2 

-y  {(y  )    +  (y  )  ) 

This  was  integrated  over  [0,  20]. 


y(o)  = 
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Table  VII.  Numerical  Results  for  Problem  5 


e 

max 
order 

steps 

func 
evals 

back 
solves 

LU 
decomps 

accurate 
digits 

time 

DIFSUB 

ID"2 

4 

50 

323 

210 

28 

-0.4 

73 

ID"3 

5 

81 

345 

264 

20 

0.4 

91 

IO"4 

5 

95 

453 

324 

32 

0.8 

115 

ID"5 

5 

140 

424 

403 

5 

1.6 

130 

ID"6 

6 

161 

487 

462 

6 

3.1 

158 

ID"7 

6 

233 

681 

656 

6 

4.2 

224 

ID"8 

6 

294 

865 

840 

6 

5.1 

283 

_Q 

10  3 

6 

401 

1174 

1149 

6 

6.1 

383 

io-10 

6 

588 

1662 

1637 

6 

7.3 

534 

ble 

nded  DIFSUB 

in'2 

5 

59 

278 

410 

18 

1.4 

100 

IO"3 

7 

63 

212 

366 

7 

2.8 

100 

IO"4 

8 

77 

247 

436 

7 

3.7 

113 

10"5 

8 

99 

307 

556 

7 

4.5 

145 

IO"6 

10 

106 

335 

596 

9 

5.1 

165 

IO"7 

10 

128 

394 

714 

9 

6.0 

228 

IO"8 

10 

157 

470 

866 

9 

6.9 

273 

IO"9 

11 

170 

505 

928 

10 

8.2 

299 

IO"10 

11 

200 

579 

1076 

10 

9.5 

361 

Again  there  is  a  significant  difference  in  performance. 
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7.      IMPLICIT  DIFFERENTIAL  EQUATIONS 

One  of  the  advantages  of  the  Nordsieck  representation  is  that  a 
predicted  value  p'   for  y'   is  available,  and  hence  implicit  differential 
equations  F(t,  y,  y')  =0  can  be  treated  without  additional   complication 
(cf.  Gear  [10]).     We  simply  chose  A     to  satisfy 

F<V  Pn  +  *0  V  P;  +  h_1  Vn1  =°-  (7J) 

The  equation  F(t,  y,  y1)  =  0  implicitly  defines  y'   as  a  function  f(t,  y) 
of  t  and  y.     Differentiating  F(t,  y,   f(t,  y))   =  0  with  respect  to  y  yields 

8F   8F   3f 


3y    3y«  9y 

and  so  the  blended  formula  has 


=  0, 


l  =  a  +  y  h  K'1  J  b  ,  (7.2) 


where  J  ~  —  and  K  " 


ay  u,,u  ,x   ay'  ' 
The  Newton  iteration  for  (7.1)  is 

A(o)  =0- 

Vl)  =  A(m)  "  [(3F/^'(m)  *0  +  h"Wy')(ll)  *irl  F(m) 

where  for  a  function  the  subscript  (m)  denotes  evaluation  at 

(t,  p  +  £~  a,   v,  p'  +  h       ju  A,    J.     The  modified  Newton  iteration  is 
v   '  K        0     (m)  1     (m) 

A/    .!\  =  A,    v   -  [K  £,   +  h  J  JLn]       h  F/   x 
(m+1)         (m)       L       1  0J  (m) 

or  equivalently 

Vn  '  4w  ■  [V  h  K"' J  V"1  K_1  hFw  • 

From  §4  we  have  a  =  [e,  1,  ...]  and  b  =  [1 ,  a,  . . .]  ,  and  hence  from  (7.2) 


we  get 


£1  +  h  K"1  J  £Q  =  1  +  (Y  a  +  6)  h  K"1  J  +  Y  (h  K"1  J)2  , 
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-1       2 
which  we  approximate  by  (1   +  c  h  K      J)    .     The  approximate  modified  Newton 


iteration  is 


Vl)  =  4W  -  [1  +  C  h  K_1  J]"2  *_1  h  r\m) 

or  equivalently 

A(m+1)   =  *(m)   -[^chJ]"1   KtK.chJ]"1   hFw    . 
Therefore  the  blended  method  can  accomodate  implicit  equations  if 
the  corrector  step  is  modified  as  follows: 

A  «-  0;  A  «-  0 

for  m  -<-  1 ,   2,   3 
r. 


F  -  F(t,  y,  y') 


r,. 


if  desirable  then  ( 


do( 


K  ■»■  3F/9y' 
cH  +■  c  x  h 
W  «-  K  +  cH  x  3F/ay 
ydecompose  W 

solve  W  x  d  =  -  h  x  F  for  d 

solve  W  x  d  =  K  x  d  for  d 


v. 


Thus  the  only  additional  work  would  be  the  computation  of  K  and  the  multi- 
plication of  a  vector  by  K.  Also  additional  storage  is  required  for  K. 

It  is  worth  noting  that  the  algorithm  does  not  depend  on  K  being 
nonsingular,  and  so  the  algorithm  should  be  suitable  for  differential- 
algebraic  systems. 
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8.  CONCLUSION 

Limited  numerical  evidence  suggests  that  the  blended  formulas  may 
be  as  good  as  the  backward  differentiation  formulas  for  stiff  problems, 
better  for  nonstiff  problems,  and  much  better  for  stiff  oscillatory  problems 
Furthermore  they  can  be  incorporated  into  existing  codes  such  as  DIFSUB 
[9],  GEAR  Rev.  3  [12],  and  EPISODE  [3],  which  have  benefited  from  con- 
siderable computational  experience.  It  is  of  interest  to  note  that  the 
AMF-BDF  blend  may  be  at  least  a  partial  solution  to  the  problem  [14,  p.  552] 
of  designing  a  method  which  automatically  selects  either  an  Adams-Moul ton 
formula  or  a  backward  differentiation  formula  for  each  individual  equation 
in  the  system. 
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APPENDIX  A. 

COEFFICIENTS 

OF  THE 

FORMULAS 

a<6> 

a<2> 

a<3> 

«W 

a<5> 

a<7> 

1 
2 

5 
12 

3 
8 

251 
720 

95 
288 

19087 
60480 

1 

1 

1 

1 

1 

1 

1 
2 

3 
4 

11 

12 

25 
24 

137 
120 

1 
6 

1 
3 

35 
72 

5 
8 

1 
24 

5 
48 

1 
120 

17 
96 

1 
40 

1 

720 
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a'<8> 

a<9> 

,(10) 

a(]1> 

,(12) 

5257 
17280 

1070017 
3628800 

25713 
89600 

26842253 
95800320 

4777223 
17418240 

1 

1 

1 

1 

1 

49 
40 

363 
280 

761 
560 

7129 
5040 

7381 
5040 

203 
270 

469 
540 

29531 
30240 

6515 
6048 

177133 
151200 

49 
192 

967 
2880 

267 
640 

4523 
9072 

84095 
145152 

7 
144 

7 
90 

1069 
9600 

19 
128 

341693 
1814400 

7 
1440 

23 
2160 

3 

160 

3013 
103680 

8591 
207360 

1 
5040 

1 
1260 

13 
6720 

5 
1344 

7513 
1209600 

1 

1 
8960 

29 
96768 

121 

40320 

193536 

1 

1   / 
72576 

11 

362880 

272160 

1 

11 

3628800 

7257600 
1 

39916800 
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(1) 


(2) 


(3) 


(4) 


(5) 


(6) 


1 

1 

1 

1 

1 

1 

1 

3 
2 

11 
6 

25 
12 

137 
60 

49 
20 

1 
2 

1 

35 
24 

15 
8 

203 
90 

1 
6 

5 
T2 

17 

24 

49 
48 

1 

1 

35 

24 

8 

1 
120 

144 

7 
240 

1 

720 


1 

363 


11 

1 
761 


44 
,(9) 

1 
7129 


140 

280 

2520 

469 

29531 

6515 

180 

10080 

2016 

967 

267 

4523 

720 

160 

2268 

7 

1069 

95 

18 

1920 

128 

23 

360 

9 
80 

3013 
17280 

1 

13 

5 

180 

960 

192 

1 

1 

29 

5040 

1120 

12096 

1 

1 

40320 

8064 
1 

362880 


(10) 


7381 
2520 

177133 
50400 

84095 
36288 

341693 
362880 

8591 
34560 

7513 
172800 

121 
24192 

11 

30240 

n 

725760 

1 
3628800 


(ID 


83711 
27720 

190553 
50400 

341747 
129600 

139381 
120960 

242537 
725760 

1903 
28800 

10831 
1209600 

n 

13440 

1 
20736 

1_ 

604800 

1_ 

39916800 
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APPENDIX  B.   SOURCE  CODE  FOR  BLENDED  DIFSUB 


IMPLICIT     REAL*8( A-H.O-Z) 

CCMMOK    /UmTGER^L    IPRPIF  .NFfN.NC  TNV.NSQLVE 

CCMMON  /INMAIN/  ICUIT 

CALL  UNDERZI  'OFF')  

£E*D ( 5.800. ENC=2C  )   I  F FOB  .  IPO i£R 
*RITE(6.lOCC)  IPROB         


1000  FOPMAT(/////'l»,<ieX,'P  R  0  E  L  E 
fcRITE (6  •100-21 


1002  FORMA T(/»0*,3X,«  TGLERANCE* «2X.»* STEPS* • 2X • • *FCN • ♦ 2X. 

1  **OECOMP*  ,2X.  **SQLVE*^2J£*  *MAX  ORDER*. 2X. 

2  'MAX  ERROR  < TC • TF ) • »3X • • T  (MAX  ERR0R)*,3X. 

3  *CPU     T  IME  •  .2X, •KFLAG* .2Xj 'STOP    AI    IF  V  ) 
IGOIT     =     0 

CO     1C     I    =    2.IFCJIER ._ 

EPS  =  10  .D0**<-I  ) 
CALL  TESTDS(EFS) 
IF  (   IQUIT.NE.O   )  GO  TO  5 
10  CONTINUE 
GO  TO  5 

20  STOP  

800  FCFMATI2I3) 
ENC 
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SUBROUTINE  TESTDS(TOL) 

IMPLICIT  ALflL».fi(  A-H.fl-/  i 

REAL*4  PW 

DI^ENSIGN  NDIM(IC) .TOJJ.O) .TFfl01.Y0<2O.lQl.v< 13f2Q>t FRRnP<?Q ). 

1  SAVE(ie.23)iYMAX(20).Pi(20.20) ,QMEGA(20 ) ,OMEGA2(20) t 

2  I  F  (  20  J  

CCKMCN  /INTGEH/   I FHCE  .  NFC N . NC I N V . NSCL VE 

COMMON  /INMAIN/ I  QUIT 

CCMMCN  /INDIF/  MAXNQ 

CATA  TO  /  O.D3.0.DQ»C«DC.O«DO.O«QO>Xl«nO  .0.00.3*0.00  / 

2CDO,1G00.00.15.DO,15.DC20.D0.2C.D0.10CO.D0.3*O.DO 

1 .00. 15*0 .DO*     - 1.00*^1 .00.-1. 00 ♦- 1^00^16 *0^D0» 

1  COO.  19*0.  DO.     2  .DO.  1  .00  .2.  CO  .17*0.  DO. 

1  .QO.C.OO.C.nC.l .00.16*0.00. **» .nn.i A*nTnn. 


CATA 
CATA 


TF 

YO 


60*0.00  / 


!  0.D0.-2.O0.-1 .00.-1 .00 .16*0. DO. 

CATA  NOIM  /   1.4,1,3.4.6.4.3*0  / 
CATA  FCURU  /  Z334 COOC OOOOOOCO 0  / 
CATA  MAXIMUM  /  1CCC  • 
NFCN  =  0 

*AXNC  =  0  

NCSTEP  -  0 
NSCLVE  =  0 
KOINV  =  0 
ERFMAX  =  O.CO 
T  MXERR  =  COO 

NECN  =  NDIM(IPROE)    __   _ 

T  =  TO(IPRCE) 

ITIME  =  0 

CC  10   I  =  1 .NECK 

Y(  1  . [  )   =  Y0< I .  IPRCB) 

YMAX(I)  =  DMAX1 { I .DC.CABSi Y<1  .1  )  )  ) 
F.MIN  =  F0URU*QMAX1(0AFS<  T)  .CABS 
F  =  TCL 
JSTART  =  0 

IF  {   T.GE.( TF(  IPRCB)  -  HMIN)   )  GO  TO 
hMAX  =  TF(IPROB)  -  T 
F  =    DMAX1 ( FM  IN.DMIN1 (HMAX.H)  ) 

NCSTEP  =  NOSTEP  -k  I 

IF  (  NCSTEP. LE.MAXNU*  )  GO  TC  35 

IQUIT  =  1 

GO     TO    50 
35     CALL     CIFSUB(NEQN.T.  Y  .SA  VE.H  .HM1  N»FWAX.TC4_^*V*1AX.  ERROR. 
1  K FLAG. JSTART.  11, Pte, IP .OMEGA .OMEGA 2) 

CALL    CALERRIT^JY-^VtiAX.GLEEBR  ) 

IF  (  GLBEWR.LE. ERRMAX  )  GO  TO  40 

ERRMAX     -     GLBERfi 

T     MXERR     =     T 
JSTART    =      1 
IF     (     KFLAG.EQ.l      ) 


10 


30 


70 


40 


GO  TO  30 

5-0  *AITE(6.9Q9)  TQL  .  NOST  EP  .  NFCN  .  NO  INV  .  NSQLVE  .  M  AXNQ  .  EflCM  AX  i  T  MX  EBP. 
1  ITIME.KFLAG.T 

999  FORMAT ( »0 • .3X.D9.3.2X.IS.2X.XC 3X.I4.SX.I4.7X.12. 
1  8X.D11.5.7X.D11.S.5X.I6.4X.I2.4X.010.4) 

RETURN  

70  IftRITE (6.999  )  TOL  .  NOSTEP . NFCN . NO INV. NSOL VE .MAXNQ .ERRMAX.  T  MXERR. 

1         IT  I*E 

RETURN 

END  ____   _. 
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SUBROUTINE  CECCMP ( P* • h • M « I P ) 
T  REAL*8< A-H.Q-Z) 


RE^L+4  PW 

CIMENSICN  PI(M)«IP(»i) 

CCMMON  /INTGER/  I PROE  .NFCN .NO  I NV .NSOLVE 
NCINV  =  NOINV  +   1 
IP(N)  =   1 
Cd  6C.K.s__l^jj 


IF  (  K.EC.N  )  GC  TC  5C 

KP1  =  K  ♦   1  - 

KM1M  =  (K  -  1 ) *M 

IPIVCT  =  K  

CO   10  I  =  KP1.N 

IF  (  ARSCPJUC  I»KM1  W)   )  .fiT  .  ARM  P*<  I  P I VHT+KMl M)  )   )  IPLVOT  s._l_ 
10     CONTINUE 

IP(K)  =  IFIVOT 

IF  (   IPIVCT. NE.K   )   IP(N)  =  -  IP(N) 

IPKM  =  IPIVCT  ♦  KM1* 

KKM  =  K  ♦  KMIM 

TEMP  =  PW(IPKN)     ' 

PW< IFKM)  =  PX( KKM) 
PW(KKM)  =  TEMP 

IF  (  TEMP.EQ.O.DO  )  GC  TO  50 
CO  20  I  =  KPl.h 

IKM  =   I   ♦  KMIM 

20        PWIIKM)   =  -  PW<  IKK|/TFMP 

CO  40  J  =  KP1,N 

J  VIM  =   (J  -  1  )*M 

IPJM  =  IPIVOT  ♦  JM1M 

KJM  =  K  ♦  JM1M  .   

TEMP  -     FW( IPJM) 

Pl(  IPJNLl.  5„EmiKJM  > 

Pto(KJM)  =  TEKP 

IF  (  TEMP. EQ.  CDC  X     GO  TH,AO^ 

CO  30  I  =  KPl.N 

IJM  =  I  ♦   JM1M         

30  PW(UM)  =     Ptt(IJM)   +  P*(  I+KM1M)  *TEMP 

40     CONTINUE 


5C     IF  (  PM(KKM)  .EC.O.C  )   IP(N)  =  0 

60  CCNTINUE  

RETURN 

ENC  
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SUBROUTINE  SOLVE (PM . N • M« IP , B ) 

U4F1_ICLT    RFAI  *fl-    <A-H.C-7) 

REAL*4  PW 

CIMENSION  P.(M),IP(«),B(M)    

CCMMON  /INTGER/  I  PROB  .NF CN .NC I NV • NSOLVE 

NSCLVE  =  NSCLVE  *  1         — 

IF  <  N.EQ.l  )  GC  TO  30 

=-_4 


CO     1C     K     =     I ,KM1 
KP1     =     K     ♦     1 
KM1M     =     (K     -     1  i  *t* 
IPIVQT     =     IP(K) 
TEMP    =     B(IPIVOT) 
EU  LP  I VOT  )     =     e  (  K  ) 


E<K)     =     TEMP 
DO     10     I     =     KP1,N 
10                 E(I)     =     B(I)     ♦    P*<  I  +  KM1M)*TEMP 
CO     2C     KB     =     1,NM1 
KM1     =     N     -     KB 
K     =     KJ«1      *     1 


KM1M    =     <K     -     1 ) *M 
B(K)     =    B(K)/PW (K+KMLM) 
TEMP     =     -     E(K) 
DO     20     I     =     1.KM1 
20  E(I)     =    6<I)     ■*■    P*(  I  +  KM1M)*TEMP 

30    Ell)     =    Edl/PHlll 

RETURN 
END 
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SUBRCUTINE     D I FFUN ( T . Y .0 Y ) 
UUJIJXULAF  „L*fl(  A-H.r>Z  ) 


CIMENSION     Yd3,2G)  ,DY(2C) 

CCMMON  /INTGEay  I  PROB  .J-tFCN.  NO  XJNV  .NSOLME 

NFCN  =  NFCN  ♦  1 

GO  TO  (10. 2C, 30. 40, 50. 6C. 70. aC. 90. 100). I PROB 
10  DY(1 )  =  -  Y( 1 ,1 i 

£E  IJJRN 


2C  ZTEMP  =  .5DC*(Y<1,1)  4  Y(l,2)  ♦  Yd. 3)  ♦  Yd, 4)) 

21  =  ZTEMP  -  Yd.lJ         

22  =  2TEMP  -  Y(  1  .2) 

23  =  2TEMP  -  Yd. 3) 

24  =  2TEMP  -  Yd.4) 

21  =  Z1*(Z1  -  LOCO. DC  ) 


22  =  22*(22  -  8CC.DC) 

23  =  23*(Z3  ♦  10. DO) 
Z4  =  24*(Z4  -  .0C1D0) 

2TEMP  =  .5DCMZ1  ♦  Z2  ♦  23  +  Z4  | 
CY(1)  -  ZTEMP  -  21 

CY(2)  =  ZTEMP  -  22 

CY<3)  =  ZTEMP  -  23 
CY(4)  =  ZTEMP  -  Z4 
RETURN 
30  CY<1)  =  -  2C0.D0*Y(1. 1  J  ♦  2000. DQ  -  (1991. DO  ♦ 
1           199.DC*T)*DE>P<-T) 
RETURN 


40  CY(1)  =  -  .1D0*Y(1,1)  -  49.9C0*Y(1.2  ) 

CY(2)  =  -     5C.DC*Y< 1 ,2) 

CY(3)  =  70.CO*Y(1,2)  -  12O.D0*Y< 1.3) 

RETURN 

50  CV(1  )  =  Y(  1  .3) 

CY(2)  =  Yd. 4)         


CY(4)  =  (  Y(  1,  1  ) *Y(  1  .  1  )  ♦  Y(  1  ,2)*Y(1  .2))**1  .5 

CY(3)  =  -  Y(  1  .1 )/DY(4) 

CY<4)  =  -  Y(l  ,2)/DY(4  ) 

RETURN 
60  CY(1)  =  -  lC.D0*Yd.l)  ♦  100.00*YC 1 .2) 

CY(2)  =  -  1  CO^J-a-YXl-d  .  -  l0.nO»YC1.2) 

CY(3)  =  -  4  ,D0*Y( 1.3) 

CY(4)  =  -  Y(l. 41  

CY(5)  =  -  ,5DC*Y< 1,5) 

CY(6)  =  -  .  lD0*Yd.6)  __. 

RETURN 
7C  ZTEMP  s  >5DC*(Y(1.1)  ±    Yd. 2)  +  Yd. 3)  ♦  Yd. 4)) 

21  =  2TEMP  -  Y( 1  ,1  ) 

22  =  ZTEMP  -  Y(i .2) 

23  =  ZTEMP  -  Y(  1  ,3) 

24  =  ZTEMP  -  Y(l. 41  

HI  =  .5D0*(Z1*ZI  -  Z2*Z2| 

_2  _-_21*-22 

Tl  =  Zl 

Zi  =  1C.DC*T1  ♦  10.OC+22  ♦  _-l 

22  =    -     1C.D0*T1  ♦  10.D0*Z2  ♦  *2 

23  =  23*123    -  1000.DC1      

Z4  =  Z4*<Z4  -  .01C0) 

ZTEMP  -  _5QX-(Z1  +  12     +  12     ♦  Z4  ) 


CV(1 )  =  ZTEMP  -  Zl 
DY(2)  =  ZTE*P  -  22 
CY(3)  =  ZTEMP  -  Z3 
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CY(4)     =     ZTEKP    -     24 

RETURN 


80  RETURN 
90  RETURN 
100  RETURN 
END 
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SUBROUTINE  CALEWR(T,Y,Y*AX, GLEERR) 

IMPLICIT  REAL*~a^A^fcUQ-Z) 

CIKENSICN  Y (  1  3.2C  )  .  Y*AX(20 ) 

COMMON    /INTGER/     IPROE •NFCN*NGJNV*NSOLV£ 

GO     TO     (  10  .  20.  30  ,  10,  50,60  ,  70,30.90.  100)  ,  IPRQB 
10    GLEERR    =     CAESdY(l.l)     -    OEXP (  -X±4^YMAXi4*4— 

RETURN 
20     Tl     =    CEJCP  C=— 1.0  0  0  .  C 0  *  T  ) 

T2 

T2 

T4 

21 

Z2 

23 

24 

ZTEMP 

GLEERR 


OEXP(-  8C0.C0*T) 
DEXP(-  1C.0C*T) 
DEXP(-  ,C01C0*T) 
10C0.DC*Tiy(Tl  -  1001*001 
800.DC*T2/( T2  -  eci.oo) 

-   io«dq/(i«-CO-*_9.co*T3) 

.C01O0*T4/( T4  -  1.C01DC) 
=  .5DC*(Z1  ♦  Z2    +  23  ♦  24) 

((Y(l.l)  -  ZTEMP  ♦  Z1)/YMAX(1))**2 
(<Y(1.2)  -  ZTEMP  ♦  Z2)/YMAX<2> )**2 
((Y(l*3)  -  ZTEMP  ♦  Z2)/YMAX(3) )**2 


(  (  Y  (  1 ,  4j  -  ZTEMP-*.. 

DSQRT(GLBERP) 


74)/YMAX(4) >♦*? 


GLEERR 

RETURN 
30  GLEERR 
1 

RETURN 
40  Tl  =  DEXP(-  5Q.CC*TJ 


DAES((Y(1,1)  -  10. CO  ♦  (1C.00  4  T)*DEXP(-T) 
1C.D0*DEXP(-  2CC.DC*!) )/YMAX(l)l 


GLeERR  = 

!  ♦ 


5C 


GLBERR 

RETURN 

CCS  T  =  DCCSdl 

SIN  T  *  DSIM T) 


((Y(l,l)  -  DEXP(-  .1C0*T)  -  T1J/YMAX(1))**2 
(<Y(1.2)  -  T  1)/YMAX(2) )**2 

((Y(l,3)  -  Tl  -  DEXP(-  120.C0*T>  )/YMAX(3  )  )**2 
DSQRT(GLSERR) 


CCS 
SIN 

SIN 
CCS 


T  J/YMAXC  1  J  J**2 
T)/YMAX(2) )**2 
T1/YMAXC3) 1**2 
T)/YMAX(4) )**2 


CLEERR    =     (  (Y  (  1.1  ) 

1  ♦     (<  Y(  1  .2) 

2  «■     { (Y(  1  ,31 

3  +     (<Y(  1  ,4) 

GLEERR    =    DSCRT(GL8ERiU 

RETURN 
60     Tl     =    DCOS( 1 C0.DC*T) 

T2    =     DSIN( 1C0.DC*T) 

T3     =     DEXP<-     13.C0*T) 

T3*(T1     ♦     T2  )  )/YMAX(  11  )**2 
T3*(T1     -     T2) )/YMAX(2) 1**2 
DEXP(-    4. C0*T)  )/YMAX(3) 1**2 
CEXP(-T) )/YMAX(4) )**2 
DEXP(-     ,5C0*T) )/YMAX(5) )**2 
DEXP(-     ♦lCO-*T))yYMAJl(6)  1**2 

GLBERR    =    OSQRT(GLBERR) 

RETURN 

70     Tl     =     OEXPC-     10.DC*T) 

T2     =    DEXP(-     1C0G»C0*T) 

T3     =     DEXP(-     ,01DO*T) 

ECCSBT    =    OCCS(-1C.OO*T) *T1 

ESINBT    =    OS  IM-1G  .D0*T1  *T1 

1  ♦     (-    9.D0*ESINBT    -     1 0  .D0*ECOS8T ) * *2 

Zl     =     -    2O.D0*(l.CO     ♦     19.DC*ECCSBT    -     ESI NBT ) /OENCM 
12     =     20.DO*(1.00     -    ECCSBT     -     1 9  .  CO *ES  INBT ) /DENCM 


GLEERR 

= 

(  ( Y(  1,1  ) 

1 

♦ 

(  (Y(  L+2-* 

2 

♦ 

(  ( Y(  1.3) 

3 

+ 

((Y(l, 4) 

4 

■f 

( ( Y(  1.5) 

5 

♦ 

( (Y(  1,61 
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23    =     1000.D0*T2/(T2    -     1001.00) 
Za     =     .0IQQ4T3/CT3    -     1.01CQ) 


ZTEMP     =     .5D0*(Zl     ♦    Z2    ♦     23    ♦     24) 

GL8ERR    =     CXYC  1  •W   ^_2TEMP     »    2  1  )  dXHAX  <  1  >  I  »*2 

1  ♦     (<Y<1,2)     -     2TEMP     ♦     Z2)/YMAX<2> >**2 

2  +     (  <  Y  C  1 ,3)    -    ZTEMP    ♦    Z3)/YMAX(3) )«*2 

3  ♦     (<Y<1«4)     -    2TEMP     ♦     24 J/YMAX (4 ) ) **2 
-GLBFRB    =    nsCBKf.LBFBB) 


RETURN 

8C     RETURN 

90     RETURN 

100     RETURN 

END 


SLEROUTINE    PEDER V ( T « Y ,P* . N ) 

IMPLICIT    REAL*a( A-h.Q-Z) 

REAL*4    PW(N,N) 

CIMENSICN    Y<13«20)- 

RETURN 

ENC  
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SUEROUT  INE  C  I  FSUE ( N • T ,Y • S AVE. H • HM IN .HM AX .EPS* MF • YMAX .ERROR, KFL AG • 

IMPLICIT  REAL*e  (A-H. Q-Z) 

C**************************+*+*+**************************************** 


c* 
c* 
c* 
c* 
c* 
c* 
c* 
c* 
c* 
c* 
c* 
c* 
c* 
c* 
c* 
c* 
c* 
c* 
c* 
c* 
c* 
c* 
c* 
c* 
c* 
c* 
c* 
c* 
c* 


THIS  SUBROUTINE  INTEGRATED  A  SEJ  OF  N  GROINARY  DIFFERENTIAL  FIRST 
ORCER  EQUATICNS  OVER  ONE  STEP  OF  LENGTH  H  AT  EACH  CALL.   H  CAN  BE 

SPECIFIED  BX.  IKE- USER  FCfi  EACH  STEP.  BUT  IT  MAY  £E  INCREASEO-OR 

DECREASED  BY  DIFSUB  WITHIN  THE  RANGE  HMIN  TO  HMAX  IN  ORDER  TO 
ACHIEVE  AS  LARGE  A  STEP  AS  POSSIBLE  WHXLE  NOT  COMMITTING  A  SINGLE 
STEP  ERROR  WHICH  IS  LARGER  THAN  EPS  IN  THE  L-2  NCfiM.  WHERE  EACH 
COMPONENT  CF  THE  ERROR  IS  DIVIDED  8Y  THE  COMPONENTS  OF  YMAX. 


THE  PROGRAM  REQUIRES  X 

CIFFUNIT, Y.CY) 

MATINVIPW  .N.M.J) 

FECERVC  T. Y.PW.M) 
THE  FIRST.  DIFFUN.  EVALUATES  THE 
VARIABLES  STORED  IN  Y(l.I)   FOR  I 
DERIVATIVES  IN  THE  ARRAY 
METFCO  FLAG  MF  IS  SET  TC  1  CR  2 


DERIVATIVES  OF  THE  DEPENDENT 
=  1  TC  N,  ANC  STORES  THE 

DOMY  IF     THE 
FOR  STIFF  METHODS.    IT  MUST  INVERT 


* 
_* 

* 
* 
* 
Jk 
* 
* 
* 
* 

* 
* 

* 


THE  N  BY  N 
SUCCESFUL. 


MATRIX  STORED  IN  THE  ARRAY  Pfe(M.M).    IF  THE  INVERSION  IS  * 
J  SHOULD  BE  SET  TO  1.  OTHERWISE  IT  SHGULD  BE  SET  TO  - 1 .    * 


PEDERV  IS  USED  CNLY  IF  MF  IS  I.  ANC  COMFUTES 
DERIVATIVES  OF  THE  DIFFERENTIAL  EQUATIONS  AS 
MF  PARAMETER.  . __ 

THE  PROGRAM  USES  DOUBLE  PRECISION  ARITHMETIC 
PCINT  VARIABLES  EXCEPT  THOSE  STARTING  WITH  P. 
SINGLE  PRECISION  TO  SAVE  TINE  ANC  SPACE. 


THE  PARTIAL 
DESCRI8E0  UNDER  THE 


FOR  ALL  FLOATING 
THE  FORMER  ARE 


THE  TEMPORARY  STORAGE  SPACE-  IS  PfinVinFO  BY-  THE  CALLER  IN-  THE 
SINGLE  PRECISION  ARRAY  P*  AND  THE  DOUBLE  PRECISION  ARRAY  SAVE. 
IEY033I  COMMENTS   DELETED  ********************************************* 
CIMENSIGN  Y(13.N)  *YMAX(N).SAVE(16,N)  «ERROR<N) . PW( N ) • CMEGAC N  )• 

1  CMEGA2CN)  .IP(N)  • 

2  A( 13)  ,E< 12)  .C(  1 1)  . GAMMA (  1  1  ),PERTST(  12.2.3) 

C»»  ♦*»***»***»***  *♦**♦♦***♦♦***,***„♦  t^t^JJL^JtJ^J-feJL^*^  ***»_»♦*»»»♦♦*♦♦»  ♦♦**♦*♦ 
C*  THE  COEFFICIENTS  IN  PERTST  ARE  USEC  IN  SELECTING  THE  STEP  AND  * 

C*  ORDER.  THEREFORE  CNLY  AECUT  CNE  PERCENT  ACCURACY  IS  NEEDEO,  * 

C************************** ******************************** ************* 
CATA  PERTST 

1  /  2. 0. 3. C. 4. 0.5. 0.6. 0,7. 0.8. 0.9. 0.10. 0.11.0.12.0. +  13.0. 

2  2.0. 1 2 ♦ C ♦ 24 . C . 37, fc^U53^  3 3 . 70 . Oa«gZ . 07 . 1 06 . 9 . 1 26 . 7.4AZ.A. 

2  168. £.191.0. 

3  3. 0. 4.  C.5. 0.6. 0.7. 0.8. 0,9. 0.10. 0.11  .0.  12.0. +13.0, +  14.0. 

4  12. C. 24. 0.37. 89, S3. 33. 70. 08.87. 97, 106. 9, 126. 7. 147. 4. 168. 8. 

4  191  .0. +213. 8. 

5  +1.0.1.0..5..1667..C4167, .008333. .001389. .0001984. 

6  +1.C. 1.0.2. 0.1.0..3158..C 74 07..01390. .0021 82. 
6        .00C2945..CC0034S2..CO0C03692, ,.0000003524  / 

A( 2)  /  -1  .DC  / 
8(1 )  /  -l.DC  • 
GAMMA  /  .C8578644D0..125D0  .. 1218908D0. .128499700. . 1087264D0. 

.0562 596100.  .  Q£ 75  A865DC  .  .  0810S623DO *.  O 759 9-87AOO *- 

.C7192936DO..C6857227D0  / 
C  /  .292893200. . 33 74973DC. »3335427CO**3A2732 9D0 . .3 169058D0 ♦ 
.299  29  7  1  DO. .286239200.  .2760 22 700. . 2677630D0 . . 2608834DC* 


CATA 
CATA 
CATA 


CATA 
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2  .2550426D0     / 

fXXXXmrKKXXPIFASF      RFMnUFKKKKKHXH«KKKKltKKHKK»KKKHIt  HKKKKXKlfyxitKltyHiOiy  lfy^| 

CCMMON     /INDIF/     MAXNQ 
CXXXXXXXXXXDOES    NOT     BELONG    hFBPXXXXUxmtHltKKXltHKititmninitxiiitKHititKirKyyyxMi^l 

IRET  =  1 

KFLAG  =  1  

IF  ( JSTART.LE.O)  GO  TC  140 

C*  BEGIN  BY  SAVING  INFCRMATIGN  FOR  POSSIBLE  RESTARTS  AND  CHANGING 

C*  H  EY  THE  FACTOR  R  IF  THE  CAILFR  HAS  CHANGE!?  H. ALL  VARIABLES 

C*  DEPENDENT  ON  H  MUST  AL SC  BE  CHANGEC. 

C*  E  IS  A  COMPARISON  FCR  ERRURS  OF-  THF  CURRENT  ORDER  NQ*  EJUS*— IS.   ._   

C*  TO  TEST  FOR  INCREASING  TFE  ORDER,  EDWN  FOR  DECREASING  THE  ORDER. 

C*  MM£M   IS   THE  STFP   S  I  7F   THAT   teAC  |  CFO  DM   THF  IA<^T   rAI  I  . 

100   D0  110I=1.N  

CO  I  10  J  -  1  .K 
110        SAVE(J.I)  =  Y(J.I) 

FCLD  =  HNEW 

PACUM  =  1.0       

IF  (h.EQ.HOLD)  GC  TO  130 
120   RACUM  =  H/HOLD  

IRET1  =  1 

GC  TO  750 
130   NQCLD  =  NQ 

TOLD  =  T         : 


IF  ( JSTART.GT.O)  GO  TC  250 
CO  TO  170 
140    IF  ( JSTART.EQ.-l )  GO  TO  160 

Q** ************************************************************** *****jl 

C*    ON     TFE    FIRST    CALL.     THE     CfiDER     IS     SET    TO     1     AND     THE     INITIAL 

C*    DERIVATIVES     ARE    CALCULATED. ! 

c** ********************************  ************************************ 


NQ  = 

1 

N3  = 

N 

M  = 

N*  16 

N2  = 

M  +  1 

h.4  = 

N**2 

N5  = 

M  «-  N 

h6     = 

N5  ♦  1 

CALL 

DIFFUM 

CO     150     I     =     l.N 
Ml     =     Nl      ♦      I 

150 YJ2.I  I  =  SAVE(  Ml  .1  )*H 

FNEW  =  H 
K  =  2 
CC  TO  100 
C***** ********  ************  ************************************** 

C*     REPEAT  LAST  STEP  BY  RESTORING  SAVED  INFORMATION. 

C****»*  *  *******«*********««***************i*****t*****ii*it*A***tt*ttitll 

160    IF  CNQ.EQ.NGOLD)  JSTART  =  1 

T  =  TCLD  J 

NQ  =  NGOLD 

K  -  NQ  «•  1  rl 

GC  TC  120 
C  * ****** ** A *********  A A ****** A ****** A** A ******  it**** A* A* *********** A ***i> 
C*  SET  THE  COEFFICIENTS  THAT  DETERMINE  THE  ORDER  ANC  THE  METHOD 
C*  TYFE.    CHECK  FCR  EXCESSIVE  CRDER.  —THE  LAST  TWO  STATEMENTS  OF 
C*  THIS  SECTICN    SET  I fcEVAL  .GT.O  IF  P*  IS  TO  BE  RE-EVALUATED 
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C*     BECAUSE     OF    THE     ORDER    CHANGE.     AND    THEN    REPEAT    THE     INTEGRATION  * 

C*     STEP     IF      IT     HAS    NOT    VET    BEEN    DONE     (IRET     =     1)     OR     SKIP     TO     A    FINAL       -  * 

C*     SCALING     BEFORE     EXIT     IF     IT     FAS    EEEN    COMPLETED     (IRET    =     2).  * 

C***** 4 **************************** ****** ************  ******************* 

170        IF     (VF.EQ.O)     GO     TO     18C 

IF     (NC.GT.ll)     GC     TO     190 

GO     TO     ( 221, 222. 222.224, 225. 226, 227,228. 229,230. 231) .NO 
IBO        IF     (NC.CT.12J     G£     TO     ISO 

GO    TO     ( 201  .202.203.  204,205,206.207,208,209.2  10.21  1,2121 ,NQ 
190        KFLAG     =     -2 

RETURN 
C******* *************************************************** ************* 

C*  TFE  FOLLOWING  COEFFICIENTS  SHOULD  BE  DEFINED  TO  THE  MAXIMUM  * 

C*  ACCURACY    PERMITTEC    EY     T  KE    KAChlhE. THEY    ARF.     list   THE  J3RDER    USED..  * 

C*  * 

C*  -1  * 

C*  -1/2.-1/2  * 

C*  -5/12. -3/4. -1/6  * 

C*  -3/8.-11/12.-1/3.-1/24  * 

C*  -251/720. -25/24. -35/72 ,-5/48 .-1/ 120  * 

C*  -95/288. -137/12C.-5/8, -17/96. -1/40.-1/720  * 

C*  - 1 9C e 7/6 C4 80. -49/40 • -203/2 70 • -4  5/  192. -7/ 1 44, -7/ 1 440  * -1/504  0  * 

C*  * 

C*  -1  * 

C*  -2/3.-1/3                                                                                                                                                               •  * 

C*  -6/1 1 .-6/11.-1/1 1                                                        * 

C*  -12/25. -7/1C. -1/5,- 1/50  * 

C*  -12C/274. -225/274. -e5/274, -15/274, -1/274  * 

C*  -ieC/441. -58/62, -15/36, -25/252, -3/252. -1/1 764  * 
C** ************************************************************ ********* 

20 1  A ( 1 )  =  -  .500 

CO  TO  2  25 

202  fi (  1  )  =  -  .4  166666666666667D0 
A(3)  a  -  .500 

GO  TO  235 

203  A( 1 )  =  -  .375C0 
A( 3)  =  -  .7500 

A (4)  =  -  .166666666666666  700    

CO  TO  235 

204  A(l)  =  -  .34861 1  1  11 1  1  1111100 
A(3)  =  -  .916666666666666700 
A(4)  =  -  .333333333322333300 
A(5>  =  -  .04166666666666667DO 

GC  TO  235  . 

205  A(l)  =  -  .329861111111111100 
A(2)  =  -  1.C41666666666667DC 
fl{4)  =  -  .486111111111111100 
A(5)  =  -  . 1 C4  1666666666667U0 
A(6)  =  -  .0C8332233332333333DC 

GO  TC  235  _ _  ..  ._. 


206  A(l)  =  -  .2  155919312169212D0 
A(3)  -  -  1 . 14 1666666666667D0 
A(4)  =  -  .62500 

A(5)  =  -  .  1770832233333333D0-  — 

A(6)  =  -  .C25D0 

a(7)  =  -  <oci38££8aaeeaeaaa9co 

GO  TO  235 

207  Ml)  =  -  .304224537037027000 
A(3)  -    -  1.225D0 
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4(4)  =  -  .751851651851851900 
JU-5-1  -=-  -  .2552083231223133D0 


Ait)     =  -  .0486111111111111100 
A(7)  =  -  . 0  C4  66 1 1 1 1 1 1 11 1111  ICO 
Ai8)     =    -  .000 19841269e41269£400 
GC  TO  235 


208        4(1)     =    -     .2948680004409171D0 
4<3).J=    -     1  .?9642fi5714?84flfinr. 


4(4)  =  -  .868518518518518500 

AC  5)  =  -  .  32  5  7£3££a£££8£aaDJI 

4(6)  =  -  .0777777777777777700 

4(7)  =  -  .0106481481481481500 

4(6)  =  -  .000793650752650793700 

4(91  =  -  .QO0024fni5fl7.ini5f  7?no 
GC  TO  235 

209   4(1)  -    -  .2869754464265714QD 


4(3)  =  -  1.3589285714285700 

4(4)     =    -     .9  76  55  423280  4233X10- 

4(5)     =    -     .417187500 

4(6)     =    -     T  l  1  l  3S A  l  fiA^fi^ft7nn 

4(7)     =    -     .01875DC 

4(8)  =  -  .00193452380952380900 
4(9)  =  -  .00011 160714285714300 
4(10)  =    -     .C0C0C275573192239£59DO 
C-C  TC  235 
210    4(1)   =  -   .? 80 189  564  4  3  936  700 


Ai3)  =  -  1.4144841269841300 

4(4)  =  -  1.C772156084656100 

4(5)  =  -  .498567C1940C35300 

4(6)  =  -  • 148437500 

4(7)  =  -  .029060570987654300 

4(8)  =  -  ,0  03  720  2  3809 52  381  DC 


4(9)  =  -  .00029968584656084700 
4(10)  =  -  .00001377865961199290-0- 


4(11)  =    -     .CCC00C275573192239€6D0 
GC  TC  235  


211        4(1)     =    -     .27426554003159900 
■Jg.^  _1  .464484  1269PA1  imp 


=  -  1.1715145502645500 
=  -  . 5793581 9QQ3S2 7 300 


.18832  26  61 55202 eDO 

.0  414303626543  21000 

.0062  11 144 79 £94 1800 
.0  006  25  2066  7  9894  18000 


4(10)  =  -  .OOC040417401528512600 
A(ll)  =  -  .COCO  01 51  56525-573  192200 


4(12)  -    -  .COC000025C521C836544170C 
GC  TO  235 


212   4(1)  =  -  .269028£467736492D0 
^    -  -4~. -53*39  186  7  ?4  38  6  700 


=  -  1.2602711640211600 
=  -  .659234  ie.2096265DQ- 


.2  3045  80  0264550300 
.05569 72461 C52  322QC 


.0094  39  48412  69  84130  0 
.001  1  192  7496  692  12  20  0 


.C0C09C939  1534  3915  34DC 
.COCA 0  4 8225 30 £64 1915  300 


.COCOOC150  3126503  1265CDO 
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A (  1  3 )  =  -  .C0C000002C676756S878681O0 

GC  TO  235  

E(  2)  =  -  l.DO 
GC  TO  201 
E<2)  =  -  1 .500 
E(  3)  =  -  .SCO 
CC     TO  202 

E(2)  =  -  I.  833333  3  3  3  3.33  333D0 

E( 3)  =  -  1 .00 

E(4)  =  -  .1666666666666667D0 

GC  TO  203 

E(2)  =  -  2.C82333233223333Q0 

E(3)  =  -  1  .4583333333233330C 

E(4)     =    -     .4  16666666d6£6i66ZQC--       

E<5)  =  -  .04166666666666667D0 

GC  TC  204 

8(2)  -  -  2.282333232223333DC 

E<3)  =  -  l.e75D0 

6(4)  =  -  .7C8333333323333300 

ECS)  =  -  • 12  500  

6(6)  =  -  .00833333333232333300 

GC  TO  205 

E(2)  =  -  2.45D0 

£(3)  =  -  2.25S555555E555560C 

E(4)  =  -  1.G20832233323333D0 

E(5)  =    -  .24305S5555555556DC 

E(6|  =  -  .C291666666666666700 

E(7)  -   -  .oci3eaEeeee£888889oo 

C-C  TC  2C6 

E(2)  =  -  2.5928571428571400 

E(3)  =  -  2.6055S5555SS5556DC 

E(4)  =  -  1  .3430S555555S556QC     - 

6(5)  =  -  .388688888888888900 

E(6)  =  -  .oe3eee£e8ee£ee889DO 

£(7)  =  -  .OC5S5555555S5555560C 

E(e)  =  -  .CCO  1984126S64126S800 

GO  TO  2C7 

E(2)  =  -  2.71785714-26571400.  „   _ 

EC3)  =  -  2.9296626984127C00 

E(4)  =  -  1.66875CC 

E(S)  =  -  .556770833332333300 

E(6)  =  -  .1  12500 

E(7)  =  -  .0  135416666666666700 

E(6)  -  -  .00C8926571 4.2B5X14 3 C 0 

E(<5)  =  -  .000024601587201567300 

C-C  TC  2C8 

E(2)     =    -     2.828968253<S6825D0 

E(3)  =  -  3.2216468253^66300 

E(4)  =  -  1  .9942680776C1410C 

E(5)  =  -  .74  2  167500  _.  _ 

E(6)  =  -  .17436342592592600 

E(7)  =  -  .0263416666666666700 

E(£)  =  -  .0023974667724667700 

£(9)  =  -  .C001240C7936507937CO 

E(10)  =  -  .COC0027557319223985900 

CO  TO  20  9  

6(2)  =  -  2.928968253S6825DC 
612)  -  -  2.5  145436507S365DO 
E(4)  =  -  2.21  7432760141C9D0 
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E<5)  =  -  .94161430776014100 
F  (  6.  >  =  -«2Afl 58 21 7592^9 2600 


E<7)  =  -  .043478009259259300 

E<8)  =  -  .00500  16E34a91EO4A0il 

E(9)  =  -  .000363756613756614C0 
E(10)  =  -  .00C 01 £156525573 192200 


E<11>  =  -  .00000027557319223965900 
24J) 


231    E(2)  =  -  3.CI 967734487735D0 

=  -  3.7808134920634900 

=  -  2.6369367283950600 

=  -  1.15  229CC13J2i-l£XDQ_ 


.3  34  1834  766  3  1393  00 
.  Q 66.026  3f  RP.t  fflPflfl90Q 
. CC895 4 1 9973544 S74D0 
.0CC818452JSaC9S23aiD0 


E(10)  =  -  .000048225306641975300 

E(ll)  =  -  .00000  165343915343915011^ 

E(12)  =  -  .C3C0000250521083854417D0 

CO  TO  £11 

235   K  =  NC+1 

ID CUB  =  K  _  

MTYP  =  (4  -  MF)/2 

ENG2  =  .5/FL0AT(NC  ♦  11 

ENC3  =  .5/FL0AT<NC  ♦  2) 

ENQ1  =    0.5/FLOAL1NOJ 

FEFSH  =  EPS 

EUP  =  <PERTST(NQ,XTYP.2)*PEPSH)**2 

E  =  (PERTST (NQ.MTYP.l  )*PEPSH>**2 

ED*N  =(PERT£T  < NG . MTYP «3 J *PEPShl±A2 

IF  (EDfcN.EQ.O)  GC  TO  780 

ENC     =     FPR*FNQ.1/riFI  nAT<NI> 

240        IHEVAL    =     MF 

GO  TO  (  250  ,  68C  )  .IRET 


C*  THIS  SECTION  COMPUTES  ThE  EBEOICIEQ  VALUES  -EUL  EFFECT  IV ELX 
C*  MULTIPLYING  THE  SAVEO  INFORMATION  BY  THE  PASCAL  TRIANGLE 

C*  MATRIX. 

Q* ************************************************************ ********* 

250   T  =  T  ♦  H 

CO  260  J  =  2.K 

DO  260  J  I  =  J.K  

J2  s  K  -  Jl  +  J  -  1 

-— CO-  26.0  L  =-  1  .N 

260  Y(J2,I)  =  Y(J2.I)  ♦  Y(J2+1.I> 

C**********4****t*********4*******4**** **********************  ********** 
C*  LP  10     3  CORRECTOR  ITERATICNS  ARE  TAKEN.    CONVERGENCE  IS  TESTED 

C*     EY  REQUIRING  CHANGES  TO  EE-LESS  THAN  END  WHICH  IVDEPENOENT  ON 
C*     THE  ERROR  TEST  CONSTANT. 

C*         THE  SUM  OF— THE-  CORRECT  ICNS  IS  ACCUMULATED  IN  THE  ARRAY 

C*  ERRCR(l).  IT  IS  EQUAL  TC  THE  K-TH  DERIVATIVE  CF  V  MULTIPLIED 
C*  BY  H**K/(FACTGfil AL(K-1  )»A<K  )  ) .  AND  IS  THEREFORS^PROPQRT IQNAL 
C*  TO  THE  ACTUAL  ERRCRS  TO  THE  LCHEST  POteER  OF  H  PRESENT.  <H**K) 
C******* ***************************************** ********m±*i 
00  270  I  =  1  ,K 

CME  6  Al  LI  -,=_  CDC 


CMEGA2I I »  =  0.C0 
270      SAVE115.I)  =  0.C0 
CO  430  L  =  1.3 
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CALL     OIFFUN     ( T . Y  •  SA VE ( N2 •  1  )  ) 

C***************  *  *  * *  *♦*  4*4  4  4**4****  **********************.*.**.*.********** 


C*     IF  THERE  HAS  BEEN  A  CHANGE  CF  CRDER  OR  THERE  HAS  BEEN  TROUBLE 
C+     IklTF  CONVERGENCE.  F*  IS  RE-EVAH*AT£D  PRIOR  TO  STARTING  THE 
C*     CCRRECTOR  ITERATICN  IN  THE  CASE  OF  STIFF  METHODS.    IWEVAL  IS 
C*     THEN  SET  TO  -1  AS  AN  INDICATOR  THAT  IT  HAS  BEEN  DONE. 

C *****  * ****************** ******** ******************************** ****** 

IF  (I  *EVA4_,LI  ♦  li-   GO  TC  350 

XMNS  HC  =  -  H*CCNQ) 
IF   (MF.EQ.2)  GO  TO  310 
CALL  PEDEFVCT. Y.PW.N3 ) 
CO  280  I  =  1.N4 
280         PMI)  =  XMNS  HC+FWCI) 

290      Ml  =  N3  +  1  __ 

M2     =    N*N11     -     N3 
CO     300     I     =     1.N12.N11 
300  P\»<  I)     =     1,0     ♦     PMC  M 

IWEVAL     =     -1 
CALL     DECCMPCPW.N.N3. IP) 

IF     CIP(N).NE.C)     GO     TO    350  _ 

GO  TO  440 
C**** ******************  4***4******************************************** 
C*  EVALUATE  THE  JACOBIAN  INTO  P to  BY  NUMERICAL  C IFFERENC ING •    R  IS  THE   * 
C*  CHANGE  NACE  TO  THE  ELEMENT  CF  Y.    IT  IS  EPS  RELATIVE  TO  Y  WITH         * 
C*  A  CIMHUN  CF  EFS**2.  * 

£***44  4*  ***************************  **  *******j**_±***  ********  ***  *********** 

310      00  32C  I  =  l.N 
320         SAVEC  14,1)  =  YC 1 ,  I) 
CO  340  J  =  l.N 

R  =  EPS*DMAX1 CEPS.CABSC SAVEC 14. J) )) 
Y(  1  .  J  )  =  YC  I  .J)  ♦  R 

C  =  XMNS  HC/R  _  _. 

CALL  CIFFUNCT.Y.SAVECN6. 1  )  ) 
DO  333  I  =  l.N 

Ml     =     I     ♦     C  J-l  )  *N3 
N12     =     N5     4     I 
N13     =     Nl     4      I 
330  PMC  Nil)     =     CSAVECM2.il     -     SAVF  C  M  .1  .  1X1*13 

340  Y(l. J)      =     SAVEC14.J) 

CO     TO     25C 
350  IF     (MF.NE.C)     GO    TC     370 

CO    36  0     I     =     l.N 
Nil     =     Nl     ♦     I 

36  0  SAVEC  14.1)     s     YC2^.I)     -    .SAVEtNiUIlAH- 

GO     TO    4  10 
370  CO     380     I     =     l.N 

Nil     =     N  5     4     I 
M2     =     Nl     ♦     I 
380  SAVECMl.l)     ■    YC2.I)     -     SA  VECN  1  2  •  1  )  *H 

CALL     SCLVE(P*.N.N3. IE.SAVE(N6. 1 ) ) 

DO    4C0     I     =     1 ,N 
N12     =     N5     +     I 
4CC  SAVEC15.I)     =     SAVE(M2,1) 

CALL     SCLVECPto.N.N3.IF.SAVEXN6.l J  I 
R     =     GAMVA CNQ )*H/XKNS     HC 

CC     403     I     =     1  *N  

M2     =     N5     4     I 

SAVEC14.I)     a    SAVECM2.1) 
4C3  SAVEC15.I)      =    R* C SA V E C  1  4  ,  I  )     -    SAVEC15.I)) 
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410         NT  =  N 
C  +  +  ********  *  *  *  *  *  *  *  *  *******i*t*ti*ttitlit*t*i**t****tt***»tti^tti*t|^^^n 

C*  CORRECT  AND  SEE  IF  ALL  CHANGES  ARE  LESS  THAN  8N0  RELATIVE  TO  YMAX. 
C*  IF  SC.  THE  CORRECTOR  IS  SAIC  TO  HAVE  CONVERGED,. 

DO  420  I  =   l.N  

Y<1. I)      =     Y<1,  I)      ♦     Ad  )*SAVEd4,  I  )     -     SAVEM5.I) 
YlPmlls.    Y(?.I)     -     SAVFC14.I)     ♦    R ( P I *SAVF C  1 5 . I  1 


OMEGA(I)     a    OMEGA(I)     ♦     SAVE(14,I) 

CMEGA2II)     =     CMEGA2CI)     +    SAVEI1S.-LI 

IF     (DABS(  SAVEd4,I)*SAVE(  15,1  )  )  .LE.  <BND*YMAX(  I>  >)     NT    =    NT    - 

420  CONTINUE  

IF     CNT.LE.0)     GC     TC     490 

430      CONTINUE     

C  +  +  **  +  *++  +  +  ***m  +  ++  +  *  +  +  +  m+*if++4*  +  +4++**  +  ++4  +  *  +  t  +  +  +  +  +  4  +  +  +  +  +  +  ++  +  m  +  **++m  + 

C*  THE  CORRECTOR  ITERATION  FAILED  TC  CONVERGE  IN  3  TRIES-   VARIOUS 
C*  PCSSIEILITIES  ARE  OECKE0  FOR.    IF  H  IS  ALREADY  HMIN  AND 
C*  THIS  IS  EITHER  ACAMS  METHOD  OR  THE  STIFF  METHOD  IN  WHICH  THE 
C*  MATRIX  PW  HAS  ALREADY  BEEN  RE-EVALUATED,  A  NO  CONVERGENCE  EXIT 

C*   IS  TAKEN.  CThERWlSE  THE  VATRIX  Pin  IS  RE-EVALUATEn  AND/OR  THF 

C*  STEP  IS  RECUCEC  TO  TRY  ANO  GET  CONVERGENCE. 

c**»*» ♦♦»♦»♦♦♦♦»♦ »♦♦♦♦» ♦»»♦♦♦ ♦♦♦♦♦♦♦♦♦♦♦*♦♦♦♦»♦♦♦ ♦»»♦♦♦♦»♦♦♦♦♦»♦♦»♦♦♦ »J 
440    T  =  T  -  H 

IF  (  (  h,LE.  (h«  IN*  1  .00CC1  )  )  .ANC  .  I  {  IWEVAL  -  M T YP  )  .LT  .- 1  )  )  GO  TO  460 
IF  (  (MF.EQ.OI  .OR. ( I WEVAL.NE  .0))  RACUM  =  RACUM*0.25D0 

IteEVAL  =  f*F  

IRET1  =  2 
GC  TC  750 
46C    KFLAG  =  -3 
470   CO  480  I  =  1,N 

CO  480  J  -  l.K 

480  Y(J.I)     =     SAVEIJ.IX 

H  =  FOLD 
NQ  =  NGCLD 
JSTAPT  =  NG 
RETURN 

C*  THE  CORRECTOR  CONVERfiFC  ANC  (TNTRni   TS  PASSED  TO  STATEMENT  520 

C*  IF  THE  ERROR  TEST  IS  O.K.,    AND  TC  540  OTHERWISE. 

C*  IF  THE  STEP  IS  O.K.   IT  IS  ACCEPTED.  IF-IODUfl  HAS  BEEN  REDUCED  _. 

C*  TO  CNE,    A  TEST  IS  MADE  TC  SEE  IF  THE  STEP  CAN  BE  INCREASED 

C*  AT  THE  CURRENT  ORDER  OR  BY  COXKG-TQ-  ONE  HIGHER  CR  ONE  LOWER-- 

C*  SUCH  A  CHANGE  IS  ONLY  MADE  IF  THE  STEP  CAN  BE  INCREASED  BY  AT 

-C-W-LEAST  -X^l^      IF  NU  CHANGE  IS  POSSIBLE  IOOUB  IS  SET  TO  10  TO 


C*  PREVENT  FUTHER  TESTING  FOR  10  STEPS. 

C*  IF  A  CHANGE  IS  POSSIBLE.  IT  IS  WACS  AND  IOQUB-J-S  SET  TO 
C*  NQ  ♦  1      TO  PREVENT  FURTHER  TESING  FOR  THAT  NUMEER  OF  STEPS. 
C*  IF  THE  ERROR  WAS  TOC  LARGE*  THE  OPTIMUM  STEP  SIZE  FOR  THTS  OR 
C*  LCWER  ORDER  IS  COMFLTED.  ANC  THE  STEP  RETRIED.    IF  IT  SHOULD 
XA-EAIL  TWT-CE-  MQRF   IT  IS  AN  INDICATION  THAT  THE  DERIVATIVES  THAT 


C*  HAVE  ACCUMULATED  IN  THE  Y  ARRAY  HAVE  ERRORS  OF  THE  WRONG  ORDER 

C*  SC  THE  FIRST  CEfllVATIVES  ARE  RECOMPUTED  AND  THE  CROFR  IS  ^SEI . 1 

C*  TC  1. 

c*»  ♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦»j.»it»»»»Aj|.jt  ♦»♦♦♦♦♦♦♦♦  ♦»»♦»»»♦♦♦  »»^t  »♦♦»♦♦»♦♦  »»jp 

450    C  at  0.0 

ERROR(I)   =  CMEGA(I)  ♦  CMEGA2(I) 

500      D  =  D  +     (ERROR  CI  )  / YMAXX  I^.U<*2 1 

IHEVAL  =  0 
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IF  (C.GT.E)  GO  TO  54C 
IF  (K.LT.3)  dO    TC  520 


C** *******  4 4************************************* *********************** 
C*  COMPLETE  THE  CORRECTION  CF  THE  HIGHER  OROER  DERIVATIVES  AFTER  A  * 
C*  SLCCESFUL  STEP.  * 

C****  *  4***4  4************************************************************ 
00  510  J  =  3.K 

CO  510  I  -  1«*     

510         Y(Jtl)  =  Y(J,I)   4  A < J)*CKEGA< I >  ♦  B  (  J )  *GMEGA2  (  I  ) 
520   KFLAG  =  +1 
FNEW  =  H 
CXXXXXXXXXXPLEASE  REMO VEX X X XX X X XX XX XXXXXXXXXX XXXXXXXXXXXXXX XXXXXXXXXXXXX 

MAXNG  =  MAXC ( MAXNC. K ) 
CXXXXX>XXXXOOES  NCT  BELONG  HERE  XXJLX  XX  X  X  X  X  X  X  XX  XXXXJLXXX-XXXXXXXXXXXXXXXXXX  X 
IF  (ICCUB.LE.l)  GO  TO  550 
IDCUE  =  ICCLE  -  1 
IF  (ID0U8.GT.1)  GO  TO  7CC 
CO  530  I  =  1 .N 
530      £AVE(16.I)  =  EPR0fi(l> 

GC  TO  700  

Q********* *******************************************************  ******* 

C*     RECOCL  THE  FAILURE  FLAG  COUNT  TO  CHECK  FOR  MULTIPLE  FAILURES.  * 

C*  RESTORE  T  TO  ITS  ORIGINAL  VALUE  AND  TRY  AGAIN  UNLESS  THERE  HAVE  * 
C*  THREE  FAILURES.  IN  THAT  CASE  THE  CERIVATIVES  ARE  ASSUMED  TO  HAVE  * 
C*  ACCUMULATED  ERRORS  SO  A  RESTART  FRGM  THE  CURRENT  VALUES  OF  Y  I S  * 
C*TRIED.  _....._..  ._.____  * 

C  +  ************* ************************************************  ********* 
540    KFLAG  =  KFLAG  -  2 

IF  (h.LE.  (HMN*1  .000C1)  )  GO  TC  740 
T  =  TOLD 

IF  (KFLAG. LE. -5)  GO  TC  72C 
C***44 4 4444 4*4***4 44* 44 4* 44 4* *********************  ********************** 
C*  PR1.  PR2.  AND  PR3   KkILL  CONTAIN  THE  AMOLNTS  BY  WHICH  THE  STEP  SIZE   * 
C*  COULD  EE  DIVIDED  AT  ORDER  ONE  LCteER.  AT  THIS  ORDER.  AND  AT  ORDER      * 
C*  ONE  HIGHER  RESPECTIVELY.  * 

C*4*4 4* ************************************************  ***************** 
550   FR2  =  (C/E)**ENQ2*1 .2 

FR3  =  1.E42C  

IF  ( (NQ.GE.WAXDER). OR. (KFLAG. LE.-l) )  GO  TO  57C 
C  =  0  .0 

CO  560  I  =  1 .N 
560      C  =  C  +     ((ERROR(I)  -  SA  VE  (  I  6  •  I  )  J  •  YMAX(  HI  **Z 
PR3  =  (D/EUP )**ENCJ*1 .4 

57C   PHI  =  1.E4-2C      -   -   

IF  (NG.LE.l )  GO  TO  59C 
C  =  CO 

CO  5E0  I  =  l.N 
580       C  =  D  +   ( f ( K. I )/YMAX{ I ) ) 44Z 
PR1  =  (C/£D*N)**ENQ1*1 .3 

590   CONTINUE  __      -  .  _..  

IF  (PR2.LE.PR3)  CO  TC  650 
IF  (PR3.LT.FR1)  GC  TC  660 
600    R  =   1  .C/AMAX1 (PR  1 , 1 .E-4) 

NEteQ  =  NQ  -  1 
610    IDCUB  =  10 

IF  (  (  KFLAG.  EQ.l  )  «AND^XR^LT-^-Ll  .  1)  )  )  GO  J-O-700- 
IF  (NE*Q.LE.NC)  GC  TC  630 
C** ********  4*** **************************************************  ******* 
C*     COMPUTE  ONE  ACCITIONAL  SCALED  DERIVATIVE  IF  ORDER  IS  INCREASED.       * 


62 

C******  ****************************************************************, 

-DC  &2.Q     La  1*M 

620      YCNE*Q*1.I>  =  EfiRCP(I)*A(K)/DFLOAT(K) 

630   K  =  NEWQ  ♦  1  

IF  (KFLAG. EC. 1)  GC  TO  670 

KACUM  =  RACUM*R       

IRET1  =  3 

_£D_IO  750 I 

640    IF  (NE*Q.EQ.NQ)  GO  TG  250 

KG  =    NEWQ  : 

GC  TC  170 
650    IF  (PR2.GT.PR1)  GO  TO  600 

NEteQ  =  NQ 

fi  =  1.0/AMAXl(ER2»l^£2=AJ 

GC  TC  610 
660   R  =  1  .0/AMAX1 (PR3.1 .E-4) 

NEteQ  =  NQ  ♦  1 

GO  TC  6  10  

670    IRET  =  2 

fi     =    DMINKR  .HMAX/r.AflF  (H>  ) 

F  =  F*R 

hNE*  =  H  r 

IF  (NO.EQ.NEtoQ)  GO  TO  68C 

NG  =  NEHQ  

GC  TO  170 
680   Rl  =  1.0 

CC  690  J  =  2.K 
£1  =  R1*R 
CO  650  I  =  l.N 
690         Y(J,I)  =  Y(J.I)*fil 

IDOUB  =  K 
700   CC  710  I  =  l^fel   


710      YMAX(I)  =  DMAX1  (YMAX(  I  )  ,DAES( Y<  1  .  I  )  )  ) 
JSTART  =  NQ 
RETURN 
720    IF   (NC.EQ.l)   GC  TC  780 

CALL  DIFFUN  ( T . Y  •  SA VE (N2 • 1 )  ) 

fi  ..a  HyFOLQ- 

CC  73C  I  =  l.N 

Yll.IJ  =  SAVEC1.I*      

Ml  =  Nl  4  I 

SAVE  (2,  I)  =  hCLD*SAV£<is.ll^U 

730      Y<2.1>  =  SAVE(2.I)*R 
Mi  _^=_4 


KFLAG  =  I 
GO  TO  170 
7AC    KFLAG  -     -1 
F.NE*  =  H 
JSTART  =  NC 
TURN 


C*****  ***************************************************************** i 

C*  THIS  SECTION  SCALES  ALL  VARIABLES  CONNECTED  WITH  H  AND  RETURNS J 

C*  TC  THE  ENTERING  SECTION.  ' 

C******************************************** ****************** ********) 

750    RACUM  =  DMAX1 (CAfcSCHM IN/HOLD) .RACUM) 

-RACUM    =     DMlKUfiACLM,CABS(HMAX/HCLO)  > 

Rl     =     1.0 


CC  760  J  =  2.K 
Rl  =  R1*RACUM 


780 
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DO  760  I  =  l«K 
760         Y(J.I)  =  SAVEC  J.I  l*&± 

h     -  HCLD*RACUM 

DC  77C  I  =  l«h 
770      Y(  1  .1  >  =     SAVEC 1 »I ) 


ICCU6 

=  K 

CZ     TC 

(   130 

KFLAG 

=  -4 

GC  TC 

470 

£NC 

25C  .  640  ),IPET1 
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